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NOMENCLATURE 


Symbol 

a 


a l 



c 


E 

F 

G 

g 

§s 


Definition 

Blade airfoil section lift curve slope, 9C^/9a 
Longitudinal rotor tilt 

2 

Body axis acceleration, m/sec 

(1) Blade flapping motion coefficient; 

(2) Covariance of innovation vector 

Lateral rotor tilt 

Fuselage aerodynamic X,Y, and Z-axis moment 
coefficients 


Hub X,Y, and Z-axis moment coefficients in 
rotor axes 


Hub X,Y, and Z-axis force coefficients in rotor 
axes 

Blade reference chord length, m (ft) 

Expected value operator 

(1) System dynamics matrix; 

(2) Statistical test function 

System noise or control matrix 

2 

Acceleration due to gravity, 9.8 m/sec (32.174 
f t/sec^ 

Structural damping coefficient 
vii 



Symbol 


Definition 


H 


w 1 * 


V 1 . 


K 

L 

L 

m 


M 

M 


M 


M. 


y 

N b 

N 

P 

P 

Q 


B 


(1) Rotor hub inplane force, along rotor X axis, 

N (lb); 

(2) System observation matrix 

^2 2 2 
Blade moment of inertia / r mdr; Kg-m (slug-ft ) 

o 

2 2 

Reference blade moment of inertia, Kg-m (slug ft ) 

Fuselage moment of inertia about X, Y, Z body 

2 2 
axes Kg-m (slug ft- ) 

Blade flap and lag moments of inertia, nondimen- 
sional 

Kalman gain matrix 

Body-axis rolling (X-axis) moment, N-m (ft-lb) 

Hub rolling moment, rotor axis, N-m (ft-lb) 

(1) Total helicopter mass, Kg (slugs) 

(2) Blade mass per unit span, Kg/m (slug/ft) 

(1) Body-axis pitching (Y-axis) moment, N-m (ft-lb) 

(2) Information matrix 

Hub pitching moment, rotor axis, N-m (ft-lb) 

Angular momentum, N-m-sec (lb-ft-sec) 

Rotor aerodynamic rolling (X-axis) moment, rotor 
axes, N-m (ft-lb) 

Rotor aerodynamic pitching (Y-axis) moment, rotor 
axes, N-m (ft-lb) 

Number of blades on the rotor 

Body-axis yawing (Z-axis) moment, N-m (ft-lb) 

State error covariance matrix 

Body-axis roll (X-axis) rate, rad/sec 

Body-axis pitch (Y-axis) rate, rad/sec 

(1) Rotor aerodynamic torque, rotor axes, N-m (ft-lb) 

(2) Process noise statistics matrix 
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Symbol 


Q 
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R 


r 


r FA 

S 

T 


V 


B 


v 



w 

X 

X B 

X HUB 



X 


I 


x 

Y 


Definition 

Rotor hub yawing moment, rotor axes, N-m (ft- lb) 
Body-axis yaw (Z-axis) rate, rad/sec 

(1) Reference rotor radius, m (ft); 

(2) Measurement noise statistics matrix; 

(3) Correlation coefficient 

Position along blade radius, m (ft) 

Blade pitch bearing radial offset, m (ft) 

Standard estimation error 

Rotor thrust (Z-axis) force, rotor axes, N (lb) 

Fuselage X-axis translational velocity, body axes, 
m/sec (ft/sec) 

Fuselage Y-axis translational velocity, body axes, 
m/sec (ft/sec) 

System measurement noise vector 

Fuselage Z-axis translational velocity, body axes, 
m/sec (ft/sec) 

System noise vector 

Longitudinal axis designation 

Body-axis X force, N (lb) 

Rotor hub longitudinal location relative to 
vehicle reference center, positive forward, m (ft) 

Rotor blade torque offset, m (ft) 

Blade chordwise bending deflection, m (ft) 

Distance of blade c.g. aft of blade elastic axis, 
m (ft) 

System state vector 

(1) Lateral axis designation; 

(2) Set of system measurements 
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System 


Definition 



y 

Z 


'HUB 


'FA 


Body-axis Y force, N (lb) 

Rotor hub lateral force, rotor axes, N (lb) 
System measurement vector 

(1) Vertical axis designation; 

(2) Blade lagging 

Body-axis Z force, N (lb) 

Rotor hub vertical location relative to vehicle 
reference center, positive downward, m (ft) 

Gimbal undersling, m (ft) 

Blade vertical bending displacement, m (ft) 


Symbol 


NOMENCLATURE (GREEK) 
Definition 


a 
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r 


Y 

6 


6 


e 




5 


F 


S TR 


6 


1 


6 


2 


6 


3 


£ 




(1) Aircraft angle of attack, rad; - 

(2) Angular displacements about rotor axes, rad 

(1) Blade flapping degree of freedom or flap 
angle, rad; 

(2) Aircraft sideslip angle, rad 

Discrete system noise matrix 

4 

Lock number, pacR /I^ 

Kroneker delta function 

Elevator, rudder, and aileron deflection, rad 
Flaperon deflection, rad 
Tail rotor collective pitch, rad 
Hub precone angle, rad 

Blade droop, outboard of pitch bearing, rad 

Blade sweep, outboard of pitch bearing, rad 
Estimate error 

Blade lagging degree of freedom or lag angle, rad 


x 






! 

Symbol 

Definition 


6 LO 

General longitudinal control deflection, rad 


6 LA 

General lateral control deflection, rad 

— 

0 

(1) Blade pitch angle, rad; 

(2) Set of parameter values 


0 B 

Vehicle pitch angle, about Y body axis, rad 


0 R 

Rotor shaft tilt relative to Z rotor axis, 
positive forward, rad 

— 

y 

Vehicle density parameter, m/pR 3 


V 

Innovation vector 


V v c 

Blade flap and lag rotating natural frequencies, 
rad/sec 

— - 

p 

Air density. Kg/m 3 (slug/ft 3 ) 

— 

a 

(1) Rotor solidity, N^c/irr; 

(2) Standard deviation 



Rotor azimuth angle, rotor axes, rad 


n 

Rotor angular rate, rad/sec 



NOMENCLATURE (SUPERSCRIPTS) 


Symbol 

Definition 


(*) 

Time derivative, d/dt 


n 

Estimated value 

. 

o , 1C , IS 

Rotor coning, longitudinal cyclic, lateral cyclic 
degree of freedom, respectively, shaft-fixed axes 


( )* 

Normalized on 1^ 

— 

( r 

Geometric derivation, d/dr 


A 

Aerodynamic 


W 

Wing 



xi 






Symbol 


Definition 


H 

V 

FUS 

TR 

T 


Symbol 


Horizontal stabilizer 
Vertical stabilizer 
Fuselage 
Tail Rotor 
Matrix transpose 

NOMENCLATURE (SUBSCRIPTS) 
Definition 


( ) 


o , 1C , IS 


b 

c 

o 


Rotor coning, longitudinal cyclic, lateral cyclic 
degree of freedom respectively, shaft-fixed axes 

Rotor blade 

Collective pitch 

Trim or steady-state condition 


xi 1 



DEFINITIONS OF INERTIAL CONSTANTS [3] 





r Z m dr/I^ 



r m dr/I^ 


Notes : 

t h 

= i l mode shape 

m = mass/unit span 

q i Mr ? 

— All terms nondimensional — 
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Cr - r FA )5 2 



I* • = 

qk<li 


L B V CZ o i_x o k - x I K) ) m dr 


f rm J n** (Z o i-x o k-Xjk) "dp dr 


/ Lg-rii (-r6 1+ (r-r FA )6 2 ) m dr 





+ / *V k B m / V* ( V' x o k ' x C krdp dr 

o o 
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(continued) 


Xlll 



m(5 FA;L " 6 FA 2 + B G^* n i} dr 

+ n k fr FA^ * ( 6 FA 7 x B + 6 FA, f 

2 3 r~ . 


-V 



• n^m dr 


UNIT CONVERSION* 


1 ft/sec 
1 m/sec 
1 "g" 

1 lb/in 2 


0.3048 m/sec 
3.2808 ft/sec 

32.2 ft/sec 2 = 9.81 m/sec 2 
703.1 Kg/m 2 6897.41 N/m 2 


^ This report supports the computer programs previously- 
installed at NASA-LaRC and NASA-Ames. These programs 
utilize English units in order to be consistent with 
the engineering units outputted by the type of data 
provided to SCI (Vt) by the government. In order to 
directly compare results with the output of these 
codes, English units have been retained herein. 
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INTRODUCTION 


The rotorcraft has presented a formidable design challenge 
to the engineering profession throughout the development of 
these unique vehicles. As the significant breakthroughs in the 
state of the art have been achieved, there has been a growing 
need for validated analytical methods in continuing to guide 
engineering decisions. Indeed, the expansion of computer-aided 
analytical methods for predicting rotorcraft aerodynamic perform- 
ance, vibrations, loads, stability, control, and handling quality 
characteristics has affected, to varying degrees, all phases of 
rotorcraft engineering. As the vehicles have become more complex, 
and the analytical methods more sophisticated, the problem 
of validating these analytical methods has itself become a 
significant issue. In particular, the validation of analyt- 
ical predictions of rotorcraft stability, control, and 
handling qualities must address a multitude of problems 
associated with vehicle test, instrumentation, and data 
processing. These validation methods must further interface 
with a wide spectrum of analytical models, from highly 
sophisticated digital simulations (e.g. C-81 and GENHEL) 
to simple transfer functions (e.g., rotorcraft handling 
quality specifications) . 

Clearly, the establishment of a systematic methodology 
for using test data to validate, and upgrade, a’ variety of 
analytical models is an evolutionary task. Significant 
advances to achieving such a methodology have been accom- 
plished for fixed-wing aircraft [refs. 1,2] and are generic- 
ally categorized as parameter identification methods. Applica- 
tions of such techniques to rotorcraft have not been extensive, 
but two notable efforts have been reported using two different 
approaches. Molusis [ref. 3] performed an extensive develop- 
ment of the Kalman filter concept to six-degree- of- freedom 
vehicle models, and reemphasized the need to include rotor 
degrees of freedom. Gould and Hindson [ref. 4] developed a 
least squares quasi- linearization procedure which recognized 
the difficulties of the sensor bias problem. These efforts 
represent significant steps toward a comprehensive rotorcraft 
model validation methodology, primarily addressing the quanti- 
fication of estimates from a given model. 

Applications of these methods, as well as results from 
fixed-wing applications, indicated that the parameter esti- 
mation solution, though an essential part of the required 



methodology, must be integrated with other requirements. 

These include: 

(1) Data quality improvement must be achieved by assuring 
the completeness and minimum accuracy level of onboard and off- 
board sensors (e.g., air data, gyros, accelerometers, force and 
moment transducers, radar, optics). Data filtering options, such 
as Fourier transform methods, must be used judiciously to avoid 
loss of desired information. Multivariable data filtering meth-; 
ods, such as Kalman filtering, should be available to calibrate 
multi- component measurements. 

(2) Size of the identifiable (validatable) model is 
inherently limited by the information content of the data and 
the computational resources available (e.g., algorithms and 
computers) . The goal is not to obtain time history matches 
of the estimated model output with the actual data (which is 
generally possible if enough parameters are included in the 
model), but rather to derive the most accurate estimates of the 
model parameters which can be obtained from the data. Even if 
a large number of parameters is identified, reasonable compu- 
tational limits must be observed. 

(3) Test design considerations , which include specific 
maneuvers to improve ' the quality of data, should be integral 
parts of an overall test objective. Test design also requires 
specification of the data processing integration to diagnose 
and resolve particular model or instrumentation problems. 

These are the principal issues which form the back- 
ground of the methodology presented in this report. The overall 
methodology has implications beyond stability and control (see 
Johnson and Gupta, [ref. 5] , but only the stability and control 
issues will be treated here. 

System identification is a generic classification for the 
methodology of integrating simulation, wind tunnel tests, and 
flight tests in order to improve or validate predictive 
models of vehicles. This classification is a significant 
extension of the process of parameter identification, which 
is but one step in the overall process, as will be introduced 
in this section. 

The following more formal definition of system identifica- 
tion has been found useful [ref. 6]: 

"System identification determines, from a 
given input/output data record of vehicle test 
response, an estimate of the physical model 
which relates the observed data." 
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This definition, implies that system identification involves 
the following: 

(1) A mathematical model : Though the governing equa- 

tions for vehicles may be known, it is frequently found that 
other equations (e.g., effects) are present, and the data 
must be analyzed to isolate what the actual model is. 

(2) Parameters of the model : This is the parameter 

identification requirement, where the coefficients of the 
model are quantified from the data. 

(3) Random errors : Random errors must be isolated and 

removed from the model effects. 

Figure l.l shows a representation of how the system identifi- 
cation procedure is used with vehicle input and output data 
to produce estimates of linear and nonlinear coeff icienf<; . 

Conceptually, the mathematical principles of system 
identification process are directly related to curve fitting 
of data to a model, and statistical analysis of the resulting 
errors. Indeed, for simple models and perfect data, elemen- 
tary least squares methods will meet most requirements for 
model quantification. The problem, particularly for rotor- 
craft applications, is that the models are not simple and the 
instrumentation is both costly and complex. The model com- 
plexity presents two difficulties, one obvious, the other 
more subtle. The first difficulty is the computer algorith- 
mic requirements imposed by processing large amounts of flight 
data and correlating these data to high order, multivariable 
models. Both feasibility and cost considerations limit the 
extent to which brute numerical force can be applied to such 
objectives. The second, and more subtle difficulty, is the 
avoidance of overparameterization of the model relative to 
the quantity and quality of data available. A frequent 
criterion of "excellence" of curve fitting techniques is the 
fit between the data and the estimated model output. In fact, 
it is now recognized that such fit can always be improved by 
adding more model parameters , even though the parameters do 
not reflect any actual physical effect (i.e., overparameteri- 
zation). While it is necessary then to obtain a reasonable 
fit to the data, it must be realized that such is not the 
explicit objective of identification. The objective of the 
identification is to achieve a validated predictive model 
which can be subsequently used for design evaluations. 




FIGURE 1.1.- INTEGRATED ROTORCRAFT SYSTEM 
IDENTIFICATION PROCEDURE 


The objective of the rotorcraft identification methodol- 
ogy discussed here is to successively treat the problems of. 
data quality, model structure estimation, and parameter esti- 
mation in a systematic fashion, iterating with the test 
maneuvers specifications to improve the identif iability of 
the stability and control coefficients from data. Each sub- 
division of the rotorcraft estimation problem has specific 
requirements for data from analysis, wind tunnel tests, or 
previous tests. The overall flow of this procedure is shown 
in Figure 1.1. 
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The key elements depicted in Figure 1.1 are integrated 
to sequentially perform the steps of: 

State estimation .— Satisfies the requirement for esti- 
mating the actual states and controls from data containing 
biases, scale factors, and random noise. 

Model structure estimation .— Satisfies the requirement 
for isolating that part of a model which can realistically 
be determined for a given set of data. The estimated model 
is the equations for the linear or nonlinear system. 

Parameter identification .— Satisfies the requirement of 
accurately estimating the parameters of the model isolated by 
the model structure estimation step. Two algorithms are 
noted, one for linear models, the other for nonlinear models, 
(This approach is taken for a computational requirements 
reduction which is realized by optimizing a program by linear 
model computations not possible with nonlinear models.) 

Input design .— Satisfies the requirement for having a 
method of specifying tests to improve the accuracy of param- 
eters and estimated models. 

In theory, all aspects of this data processing could be 
performed with one single algorithmic method. This "procedural 
ambiguity" is admittedly present. The resolution shown in 
Figure 1.1 incorporates the desire to utilize various types and 
qualities of a priori data to optimize a particular step within 
minimum cost/maximum accuracy constraints. 


Organization of Report 

Chapter 2 summarizes the rotorcraft mathematical models 
which are the basic identification models of the algorithms 
discussed in the following chapters. Chapter 3 presents the 
background, requirements and selected state estimation algo- 
rithms. Chapter 4 discusses the general model structure 
estimation problem, and details the algorithm for rotorcraft 
estimation. Chapter 5 reviews the parameter identification 
algorithms and Chapter 6, the input design requirements and 
algorithms. Chapter 7 reviews specific numerical results 
from each step of the algorithm. Chapter 8 presents conclu- 
sions of this effort. 
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CHAPTER II 


ROTORCRAFT MATHEMATICAL MODELS FOR SYSTEM IDENTIFICATION 
Specification of Rotorcraft Model Criteria 

Overview . - This section considers the analytical models 
selected for implementing the rotorcraft parameter identification 
methods employed in this study. As is well known, detailed 
dynamic models of rotorcraft tend to be of great complexity due 
to the large number of degrees of freedom and time-varying non- 
linear characteristics of parameters that must be considered. 
Given that the principal objective of this work is the identi- 
fication of parameters relating to fuselage and rotor dynamics, 
and not a study of the modeling task per se , the basic model 
requirements were fixed by questions of system identif iability 
rather than model generality. 

The principal guidelines for model selection are summarized 
in the following criteria. 

Math model criteria . - The objectives of the parameter iden- 
tif icaTTon - 17oFir^Tar[ne^ — in this effort governed the determination 
of applicability on the part of several alternative model formu- 
lations. It was established that model selection should be made 
relative to the following criteria: 

(1) The model must contain all degrees of freedom of 

/ importance to rotorcraft handling qualities analyses. 

Rotor modes must be explicitly shown and not reduced 
to quasi-static form by inclusion in fuselage equations 
of motion. The rotor modes required are the first 
harmonics of longitudinal and lateral flapping and 
lagging blade motion, plus coning. 

(2) Explicit representations must be given to both fuse- 
lage and rotor aerodynamics. Fuselage effects must 
be shown as component build-up so that component 
effects may be determined if suitable instrumentation 
is present. 

(3) The model must be usable with minor modifications for 
both linear and nonlinear analyses as would be 
required in this effort, in order to provide the con- 
sistency needed to isolate the effects of particular 
system nonlinearities. 

(4) The model must be general enough to treat different 
types of rotors and different numbers of blades. 

(5) The model should relate to a known theoretical treat- 
ment of the rotorcraft dynamics problem, in order to 
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give insight into the origins and assumptions implicit 
in the parameters identified. That is, a general- 
purpose identification model with arbitrary parameters 
should not be used; state selection should be in 
accordance with previous analytical work. 

(6) The model should utilize both dimensional and non- 

dimensional parameters in a manner that provides the 
simplest and most meaningful representation of rotor- 
craft inertial and aerodynamic characteristics from 
the point of view of the handling qualities analyst. 

The characteristics of various methods of rotorcraft model- 
ing are examined in light of these requirements in the following 
subsection. 


Math Model Selection 

Math model types .- There exist three basic types of rotor- 
craft models that have found wide applicability in analytical 
studies and have been used by many researchers; In order of 
generally increasing complexity, they are: 

(1) Quasi-static, 6 DOF models in which rotor dynamic 
effects are assumed negligible (by setting rotor state 
rates equal to zero) and rotor aerodynamic performance 
effects are combined with fuselage effects. Rotor 
forces and moments may be represented in terms of 

the classical rotor parameters inflow ratio (X) , 
advance ratio (q) , and blade pitch (9). The equations 
may be linearized and analyzed in a manner analogous 
to fixed-wing aircraft methods [ref. 7']. 

(2) Rotor blade dynamics models, in which the rotor degrees 
of freedom are represented explicitly and rotor forces 
and moments are coupled to the fuselage through con- 
straints at the rotor hub, and in which rotor blade 
aerodynamic characteristics are determined using one 

of several azimuthal averaging techniques [ref. 8 ]. 

(3) Rotor blade dynamic models such as in (2) , except 

that rotor blade aerodynamic characteristics are deter- 
mined exactly at a large number of azimuthal stations 
and accelerations and rates are integrated numerically 
from one station to the next to determine blade motion 
in time-history form. This method, implemented on a 
digital computer, is powerful in terms of high accur- 
acy and nonlinear capability. No averaging assumptions 
are required, and the accuracy is limited only by 
integration errors and computer truncation errors. 
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The first model type is very useful for vehicle handling 
qualities studies in which blade dynamic modes may not be of 
concern, as it is of low order (8 states) and contains easily- 
determined parameters. It is useful either in a nonlinear, time- 
history output form or in a linearized, small-perturbation, 
eigenvalue and eigenvector analysis form. However, the absence 
of explicit rotor blade dynamics which leads to this simplicity 
makes the model inappropriate for studying blade aerodynamics 
effects, which is one of the principal objectives of modern 
rotorcraft research. 

The third type of model contains a very high level of de- 
tail in its time history calculations, and represents in fact a 
simulation of rotor dynamics that may include as many parameter 
nonlinearities as desired due to the flexibility of digital com- 
puters in integrating the equations of motion with variable 
coefficients. The capability to perform modal analyses is lost 
with this type of model, however, due to its numerical nature. 

Modal analyses require a model expressed in linear mathematical 
form. The maximum of analytical flexibility comes with a linear 
model der ived from the nonlinear s imulation, at a given flight cond 
tion. In this linear model, the coefficients of the differential 
equations represent averages of nonlinear parameters around the 
rotor azimuth. 

This comprises the second type of model. Several efforts 
have led to analytically-determined, linear models that show the 
origins of the nonlinear coefficients that must be arranged, 
and then impose averaging techniques to arrive at a linear, 
closed-form system of equations. Outstanding among these ap- 
proaches is that of Johnson [ref. 8], 

In the case of linear parameter identification, a linear, 
state-space model is employed in the identification program. For 
identification of nonlinear parameters, the same basic linear 
model is used but additional, nonlinear effects are accounted 
for through the introduction of functions to represent equation 
coefficients and the inclusion of nonlinear terms in the forcing 
function and inertia models. The basic method of determining 
the parameters is unchanged, though the numerical operations are 
more involved. This flexibility suggested the selection of a 
model that was linear at given flight points, but that could 
be modified by additive constants and nonlinear coefficient 
functionals for successive flight conditions or nonlinear para- 
meter estimation. 

The type two math model selected is based on the general 
work of Johnson [ref. 8] specialized to the single-rotor heli- 
copter case.- This model, derived by rigorous aeroelastic analy- 
sis, contains all of the required characteristics for the present 
purpose and provides a sound framework to serve as a basis for 
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both linear and nonlinear parameter identification. The follow- 
ing subsection reviews the background of this model and its 
application to the present parameter identification effort. 


Linear and Nonlinear Math Model Description 

Overview . - This section contains a review of the formulation 
of the set of equations used as the basis for parameter identifi- 
cation work in this study. Many of the features of the model 
developed in ref. 8 were neglected as not appropriate to experi- 
mental handling qualities studies, particularly the provision 
for higher-harmonic blade bending mode shapes. Rotor dynamics 
are determined from rotor inertial and aerodynamic characteris- 
tics, transformed to hub forces and moments, and added to fuselage 
inertial and aerodynamic forces and moments, from which acceler- 
ations of and about the vehicle center of mass are then computed. 

A schematic of the flow of computations in this process is shown 
in Figure 2.1. In the following, the particular elements shown 
in this figure are discussed. 

Formulation of equations of motion . - 

(1) Coordinate Systems and Transformations . Rotor blade 
inertial and aerodynamic characteristics are calculated in a 
coordinate system fixed with the blade (i.e., rotating with rotor 
angular velocity), then Fourier-transformed into a non-rotating, 
shaft-fixed coordinate system. This transformation is well- 
discussed in ref. 9. The rotor may then be viewed as a cone 
whose semivertex angle and inclination relative to the shaft 
change in accordance with blade motions, transformed. Vehicle 
reference axes are fixed in the rotorcraft fuselage, and the 
relationship between rotor and vehicle degrees of freedom, 
forces, and moments is illustrated in Figure 2.2. Transforma- 
tions between the two orthogonal, right-handed coordinate systems 
are shown in Figures 2.3 and 2.4, and are summarized in equation 
form in Figure 2.5. In the equations used in this study, and 
presented below, axis transformations were expanded and written 
out in equation form. 

(2) Nondimens ionalizat ion . Particular care was given to 
the use of dimensional and nondimens ional parameters in the 
equations of motion. The generality and conciseness of using 
nondimensional coefficients introduces computational problems 
if points of singularity (e.g., zero airspeed in hover) occur 
and when operating conditions depart substantially from reference 
conditions. Dimensional parameters, on the other hand, relate 
more directly to flight measurement data but are scale- and 
flight condition-dependent. In this study, both types of para- 
meters were used. As will be seen below, the basic nonlinear 
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FIGURE 2 . 1 .- SCHEMATIC OF EQUATIONS OF MOTION 
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equations are dimensional, with physical units of force and 
moment, though fuselage aerodynamic contributions are expres- 
sed as nondimens ional coefficients multiplied by the refer- 
ence nondimensionalizing terms. All nondimensionalizing was 

2 4 2 5 

done on rotor characteristics: pfi R tt for forces, pQ R it for 

moments. Conversion from conventional, fixed-wing coefficient 
definitions for components (such as horizontal and vertical 
tail) must be made where required. 

The linear equations are nondimens ional , on rotor char- 
acteristics. The nonlinear equations are dimensional but 
feature nondimens ional rotor equations which are dimensional- 
ized prior to the parameter identification step. This ap- 
proach was taken to avoid complete dimensionalization of the 
complex inertial and aerodynamic constants in the basic rotor 
equations. An added factor in the equations of reference 
is that they are also normalized on rotor blade inertia; this 
was not retained in the present equations, resulting in addi- 
tional scaling changes between the two sets. 

A summary of the final nondimens ionalization definitions 
is given in Table 2.1. 

(3) Fuselage Aerodynamics . Fuselage aerodynamics are 
expressed in coefficient form based on rotor characteristics. 

A component buildup method is used, in which the force and 
moment contribution of each component ( in the presence of 
all the others ) is added to yield net total fuselage forces 
and moments . A typical buildup is shown in Table 2.2. It is 
important to note that the contributions of individual compo- 
nents cannot be separated unless the appropriate load-cell 
instrumentation is installed on the components; otherwise, 
the sum of all components will be identified as a parameter 
varying with, say, w, another parameter varying with q, and 

so on. This is shown by the arrangement of terms in Table 2.2 
and the equations below. Groups of these terms are parameters 
to be identified. 

Note that flap and control deflections are included, as 
are tail rotor force and moment. And, in the complete moment 
equations (below), the effects of rotation of the tail rotor 
(and engine) angular momentum vectors are included. 

(4) Rotor Equations of Motion . The rotor equations of 
motion describe , in shaft axes, the motion of the rotor disk 
in response to aerodynamic and inertial forces and moments. 

The rotor forces and moments determined by these equations 
are transferred to the fuselage as hub forces and moments, 
and are discussed in the following paragraphs. The rotor 



TABLE 2.1. - NONDIMENS IONALI ZING FACTORS FOR 
EQUATIONS OF MOTION 


Parameter Units Factor 


1 . 

Length 

m (ft) 

R 

2. 

Linear velocity 

m/sec (ft/sec) 

S2R 

3. 

Angular 

velocity 

rad/sec 

SI 

4. 

Mass 

Kg (slugs) 

PR 3 

5. 

Force 

newton (lb) 

PS1 2 R 4 tt 

6. 

Moment 

nt-m ( ft- lb) 

pft 2 R 5 7r 

7. 

Moment of 
inertia 

kg-m 2 (slug-ft 2 ) 

PR^TT 

8. 

Angular 

Momentum 

kg-m /sec (slug- 

pfiR^TT 



ft 2 /sec) 




(Force) 

(Moment) 

(Moment of inertia) 
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equations used in this study, derived from reference 2 and 
retaining only first-order blade flap and lag modes, are 
shown in Table 2.3. It will be- noted that single parameters, 
either constant or nonlinear functionals, are used in these 
equations to represent the complex, theoretically derived 
coefficients of reference 2. By referring to these theoret- 
ical expressions, the content and derivation of identified 
parameters may be determined. It will also be noted that 
fuselage rates and accelerations have been transformed into 
the rotor (shaft) axes. 

(5) Rotor Hub Forces and Moments . Rotor dynamics are 
determined completely by the rotor equations discussed above. 

As the result of these motions, forces and moments are im- 
posed on the rotor hub; and in this way, rotor forces and 
moments are both aerodynamic and inertial in origin. Table 
2.4 shows the equations for the hub forces and moments, in 
rotor axes. The aerodynamic parameters in these equations 
must be identified, along with those relating to rotor and 
fuselage aerodynamics. 

Time-varying rotor forces and moments have been modeled 
in this study using the averaging method of Johnson [ref. 8] , 
which replaces periodic rotor force and moment coefficients 
with constant, Fourier- averaged values of the coefficients to 
obtain a set of constant-coefficient ordinary differential 
equations amenable to rigorous mathematical analysis. All of 
the effects of the rotor on the airframe have been modeled, in 
the present study, by computing rotor-related hub forces and 
moments under this constant-coefficient approximation and 
resolving them into body axes for the usual force and moment 
summation . 

General, nonlinear vehicle equations of motion .— Using 
the coordinate system definitions, nondimensionalization 
technique, and analyses of vehicle component aerodynamics 
and rotor dynamics, a general set of complete vehicle equa- 
tions of motion may be written. The basis for the general 
equations is Euler's equations referenced to a set of body- 
fixed axes; these equations are shown in Table 2.5. They are 
shown here as dimensional equations, and the degrees of free- 
dom appear as their true values, not as perturbation quantities. 

The nonlinear vehicle equations for use in this study 
are shown in Table 2.6. They are obtained by introducing the 
expressions developed in the previous subsection for fuselage 
aerodynamics and hub reactions due to the rotor. While the 
rotor equations given in the preceding subsection are for 
small-perturbation motion from trim, their use to compute 
actual hub forces and moments is valid, as biases due to 
steady state (trim) conditions are included in all of the 
equations and are identified along with the aerodynamic 
parameters . 
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Rotor Dynamics (Continued) 
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TABLE 2.5 - GENERAL NONLINEAR EQUATIONS OF MOTION 


General Nonlinear Equations of Motion 
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This model is employed for nonlinear simulation and 
parameter identification work, for which time histories are 
the principal output required. 

Linear vehicle equations of motion .— The linear equations 
of motion for the rotorcraft, shown in Table 2.7, are a special 
set which are derived from the nonlinear equations by the 
well known linearization technique using a first-term Taylor 
expansion of all parameters and states. Note that the effects 
of gravity and controls are also linearized. This model is 
used for linear parameter identification and general model 
analysis studies, and adds the last of the modeling capabili- 
ties desired in this study. 

Modeling the Two-Bladed Rotor .— The modeling aspects of a 
two-bladed rotor differ from those of a rotor with three or more 
blades because the two-bladed rotor is not axisymmetric and its 
inertial properties, as well as its aerodynamics, possess period- 
ically-varying characteristics. The modeling of this rotor 
is again based on the method employed by Johnson [ref. 8] , in 
his analysis of the general rotor dynamics problem. This method 
substitutes constant, average values of the coefficients for 
this time periodically-varying values, and in so doing, achieves 
simplicity and clarity in the dynamic formulation. As noted 
in reference 8, it introduces some errors because the periodic 
and nonaxisymmetric properties of the testing rotor are pro- 
nounced. Nevertheless, the method, appropriate to rotor 
dynamics studies, is certainly also appropriate to rotorcraft 
handling qualities studies, and provides for continuity of 
method in the present model formulation. The constant coeffi- 
cient method is particularly well suited to the sophisticated 
system identification methods employed in the present study 
because periodic effects will be distributed throughout the 
identified model in such a way that maximum dynamic fidelity 
is obtained. In this way, the results of system identifica- 
tion analysis of the two-bladed rotor will guide future 
modeling effort toward useful model structures for handling 
qualities and flight dynamics studies. 

Identification and State Estimation in the Two-Bladed 
Rotor~ For the purpose of system identification, a linear 
time- invariant model of the two-bladed rotor gives a high 
level of modeling error at 2 per rev. This modeling error 
occurs because the two-bladed rotor possesses a plane of 
symmetry which rotates with the blades and crosses a refer- 
ence point twice in each rotation. In a steady state, this 
error term can be filtered out. The filtering is most 
effective when low frequency characteristics are desired. 
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Thus, accurate estimates of stability and control parameters 
may be obtained from models given previously. 


The problem is made more complicated because of the 
changing level of the 2 per rev error in a transient. A con- 
tinuously tracking filter is required to minimize the error. 
This is particularly significant for system identification 
where most information about parameter values is available in 
a transient. The computer programs based on the models pre- 
sented in this chapter may incorporate such tracking filters. 


When rotor-blade measurements are available, difficulties 
are caused in transformations from the rotating system to the 
fixed system. In three or more bladed rotors, the positions 
of the blade tips at any time define the instantaneous posi- 
tion of the tip-path-plane (TPP). For a two-bladed rotor, 
two points on the blade cannot define an instantaneous plane. 
The following procedure may be used, however, for the trans- 
formation. 
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Differentiating the above equation once and then twice, we get 
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since the dynamic system of equations defines the relationship 
between the tip-path-plane acceleration and the corresponding 
velocity and position, we have 
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We have six equations and six unknowns $ , 3^ , 3 , 

• • 

3ic> 6^ s - We can determine these unknowns in terms of 3^, 

3 2 , and their first- and second-order derivatives and rotor 

inputs. The transformation is both time varying and param- 
eter dependent. In addition, both flapping rate and flapping 
acceleration measurements are required. 

The problem may be simplified somewhat by taking the 
sum and difference of 3^ and 8 2 before time derivatives 

are derived. Equations for 8 lc and B ls are still complex. 

A more direct approach is based on using the blade- 
flapping measurements without transformations to the fixed 
axis system. Periodic gains with 1 per rev periodicity are 
then obtained. The implementation is more difficult because 
of the need to store a time history of Kalman filter gains. 
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CHAPTER III 


ROTORCRAFT STATE ESTIMATION 


State estimation is the process of quantifying the time 
histories of significant vehicle response variables, based 
on measurements which themselves are susceptible to errors. 
The input data to the state estimation process are the meas- 
urements and the desired outputs are best estimates of the 
vehicle states and of the measurement errors (e.g., biases, 
scale factors, etc.). 

The first section reviews the specific requirements for 
rotorcraft state estimation. The second section then summar- 
izes the fundamentals of state estimation, followed by the 
third section which presents the basis of a rotorcraft state 
estimation algorithm. The fourth section discusses the 
structure of the computer program DEKFIS (for Discrete 
Extended Kalman Filter and Smoother. 

Requirements for Rotorcraft State Estimation 

Rotorcraft flight testing for handling qualities, sta- 
bility and control, and simulation validation objectives re- 
quires the measurement of a multitude of variables. Examples 
of such measurements are on-board variables such as accelera- 
tions, velocity, position, attitude, angular rates, and 
control inputs. Off-board measurements may also be required. 
Radar tracking is an example. These measurements are, in 
general, subject to instrumentation errors due to random 
noise, bias, scale factor, channel cross-talk, and dynamic 
lags. Further errors are introduced with instrumentation 
data conditioning and transfer due to sampling effects and 
word length limitations. The determination of the best esti- 
mate of required vehicle states from these measurements is 
therefore non- trivial. 

The systematic methodology for estimating the vehicle 
states and for estimating the errors in the measurements is 
denoted as state estimation. This methodology has evolved 
very rapidly in the past 15 years as digital computers have 
become more widespread as a flight test data processing tool. 
Parallel to the use of these digital computers has been the 
development of algorithms for minimizing the error of the 
estimated time histories relative to the actual responses. 

The number of measurements required for rotorcraft dynamic 
state description presents a formidable state estimation 
requirement for these algorithms. 

The general requirements for rotorcraft state estimation 
are to provide the following improvements over raw flight 
data : 



(1) data consistency for related attitude, rate and 
acceleration measurements (e.g., between radar and 
integrated linear accelerometer measurement posi- 
tion) ; 

(2) the capability to handle redundant instrumentation 
at different locations; 

(3) corrections for instrumentation errors such as 
scale factor and bias errors, instrumentation dy- 
namic lag, cross-talk and random noise; 

(4) estimates for gusts and other unmeasured states; 
and 

(5) estimates for all states. 

The specific requirements for rotorcraft state estima- 
tion are as follows: 

(1) estimation of fuselage states in presence of gusts, 
and estimates of gusts and gust statistics; 

(2) estimation of rotor states; and 

(3) estimates of component forces and moments. 

These general and specific requirements pose a significant 
computational task because of the large number of vehicle, 
instrumentation, and statistical parameters which must be 
estimated . 

Resolution of these requirements involves three princi- 
pal issues. These issues are: 

(1) the state estimation algorithm; 

(2) the specific models to be estimated; and 

(3) the algorithm implementation. 

These issues are addressed in the following subsections. 


Review of State Estimation Algorithms 

There is now an extensive body of experience and docu- 
mentation which provides extensive background to the science 
and art of modern applied state estimation [10]. This sec- 
tion is a summary of the principal aspects of this technology 
which are necessary for understanding the rotorcraft estima- 
tion approaches discussed in subsequent sections. 
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Linear continuous and discrete system models . - The gen- 
eral problem of state estimation is based on a formulation of 
the linearized system dynamics and measurement equations, viz 


x = Fx + Gw , x(o) given (3.1) 

y = Hx + v (3.2) 

where 

x A system state 

w A system noise (random, white process with 
power spectral density, Q) 

y A^ system measurements 

v A system measurement noise (random, white pro- 
cess with power spectral density, R) 

F A^ system dynamics matrix (assumed constant) 

G A system noise matrix (assumed constant) 

H A system observation matrix (assumed constant) 

(Explicit vectors and matrices for rotorcraft application 
are given in the next section.) 

Equations (3.1) and (3.2) are a continuous model for the 
system. The corresponding discrete model, describing system 
characteristics at discrete times i and i+1 (0 A i A N) , is 

x k+l = $ x k + Fw k » x o § iven ( 3 - 3 ) 


^k = H k x k 


+ v. 


(3.4) 


where x^.. 


y^., w k anc * v k are t ^ e sam P^ e< ^ values at time 


k of the corresponding terms in Eqs. (3.1) and (3.2), and the 
matrices <J> , r, and H are the discrete equivalents of 
F, G, and H. 


State estimation algorithms .- The state estimation 
problem is the determination of - the continuous^estimate of 
x, x, or the sequence of estimates of x^, x^, from 

measurements y or yj,, respectively. The estimate x is 

obtained from Eqs. (3.1) and (3.2) as 
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x = Fx + PH T R _1 (y-Hx) , x(0) = x given (3.5) 

P = FP + PF T - PH T R -1 HP + rQT T , P(0) assumed (3.6) 
and the estimate x from Eqs. (3.3) and (3.4) as 

x k " * x k-i 

X Q - x 0 , given (3.7) 

p k = ((^Pk-i^ 1 + r Q rT 3 ’ 1 + hJr’ 1 ^]" 1 , 

P(0) assumed (3.8) 


The following characteristics are apparent from these 
equations : 

(1) The filter is initialized with estimates of the 
initial state x(0) and state error covariance, 
P(0). In general, these initializations must be 
assumed. 

(2) • Estimates of the process and measurement noise 

statistics, Q and R, must be provided to the 
algorithm. These statistics must be estimated 
a priori (or at least by a modification of the 
algorithms described by Eqs. (3 . 1) - (3 . 4) ) . 

T - 1 

The quantity PH R is denoted as the Kalman 
gain, K, 

K = PH T R _1 (3.9) 

and is explicitly dependent on R and on Q 
(Eqs. (3.6) ,(3.8)) . 

(3) The following special cases may be identified from 
the discrete algorithm equations, (3.7) and (3.8). 

(a) <t> known, Q is finite 

This corresponds to the case shown in Eqs. 

(3.7) and (3.8). 
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(b) 4> is unity, Q is finite 


Equations (7) and (8) become 



(3.10) 

p k ’ [(P k . 1 + rQr T )’ 1 +H T R“ 1 H]' 1 

(3.11) 


which represents a fitting process for a 
process where the value of the state at time, 

K, x k , is related to the value at time k-1, 

x k-l’ within the standard deviation associ- 
ated with Q. Note that this implies no 
knowledge of the deterministic dependence of 
on As Q increases, the infor- 

mation provided to x k by knowledge of x k 

decreases. This type of representation is 
denoted as a random walk model. 

(c) $ is unity, /Q approaches infinity (e.g., large) 

Equations (3.10) and (3.11) become 


x k = (H T R" 1 H)' 1 H T R‘ 1 y k 


(3.12) 


P k = (H T R" 1 H)‘ 1 


(3.13) 


which corresponds to the well known least 
square single point estimate of x k given 

the measurement y k> and is independent of 

estimates etc. 

These considerations are relevant to selection of the a pri- 
ori assumptions on Q. Similar conclusions are appropriate 
for the continuous equations (3. 3) -(3. 4). 


Smoothing . - The estimation algorithms of the last sub- 
section yield an estimate, at any time t, of the state, 
x(t), based on measurements y(t), up to and including 
time t. One would expect, however, that the accuracy of 
such point estimates could be increased if, in addition to 
the past data up to t, future data from t to T could 
also be incorporated. In the context of sequential estima- 
tion approaches (as opposed to a batch estimation wherein all 
data are simultaneously processed), smoothing is the generic 
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operation for using future data to improve intermediate time 
estimate accuracy. The Kalman filter estimation algorithm 
can achieve this by processing the data forward until the 
last data point. The filter and smoother estimates are then 
identical because they have the same information available. 
Reverse smoothing then consists of stepping back from the 
last measurement and forming the smoothed estimate by adding 
a correction term to the forward filter estimates. 

The process of stepping back through the data can be 
shown [10] to not only reduce the uncertainty in the state 
estimate, but also provide an estimate of the uncertainty in 
the state (e.g., process noise). These characteristics then 
lead to specific modes of smoother operation denoted as fixed 
point smoothing and fixed interval smoothing. Fixed point 
smoothing is the mode in which the accuracy of a specific 
point estimate is improved. In particular, the initial con- 
dition estimate of state and its error covariance are im- 
proved based on all the data. Fixed interval smoothing is 
the use of data over an interval to provide improvements to 
the estimates of the measurement noise statistics and the 
process noise statistics. 

It should be noted that smoothing type of iterations 
can be performed over any span of data. Hence, local, iter- 
ated smoothing is frequently used for improving estimates of 
states whose actual transition dynamics are nonlinear. 

The smoothing operation on data in a nontrivial com- 
putational problem; the filter requires the storage 
of estimates and covariances of, at most, two successive 
data points. The smoothing solution requires that these 
estimates and covariances be retained for all data points . 
Further, the backward propagation of the filter may, for a 
system whose transition matrix is stable, become unstable 
(depending on data rate) . Such a combination of high com- 
puter storage requirements and possible numerical instabili- 
ties has resulted in few examples of useful filter/ smoother 
implementations . 


Rotorcraft State Estimation Algorithm 

The high dimensionality of rotorcraft state estimation 
requires a highly systematic integration of the filter/ 
smoother characteristics discussed in the last section. The 
DEKFIS program utilizes a Friedland-Duffy extended Kalman 
filter with a locally iterated smoothing algorithm to pro- 
vide the filter estimates. An additional fixed interval 
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smoothing algorithm can be used to estimate the gusts (i.e., 
process noise) and provide smoothed estimates. This section 
details the specific organization of the rotorcraft state 
estimation algorithms. 

Structure of model and filter/smoother . - The fundamen- 
tal rotorcraft sensor model consists of two types of elements 
in the filter state vector. These are the time-varying 
states of the vehicle dynamic model and "constant states" 
associated with sensor biases. (These constant states will 
therefore be referenced as biases.) The general filter 
model will therefore be of the form 

x-^Ct) = f (x 1 ,x 2 ,u,t) + T(t)w(t) , X]_(0) = x 1q (3.14) 
x 2 (t) = 0 , x 2 (0) = b Q (3.15) 

with measurements 

y (t) = h(x 1 ,x 2 ,u,t) + v(t) (3.16) 

where X]. is the dynamic state initial condition and b Q 
is the a priori estimate of bias, and noise statistics 


E{w(t) } = w(t) 

E{ [w(t)-w(t)] [w(t)-w(t)] T } = Q (t) 6 (t-r) 
E{ v (t) } = v(t) 

E( [v (t) - v (t) ] [v(t)-v(t)] T > = R(t) 6 (t-r) 


Equations (3.14) and (3.16) are nonlinear and must be lin- 
earized in order to satisfy the linearity requirement of the 
Kalman filter solution. These linearized equations are of 
the form 


x^ = Fx^ + Bx 2 + Gu + Tw (3.17) 

x 2 = 0 (3.18) 

(3.19) 


y = Hx-^ + Cx 2 + Du + v 
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Using the Friedland-Duffy approach, the filter equa- 
tions can be written as sets of bias-free filter equations 
for the primary states and bias filter equations for the 
bias states. 


Bias-Free Filter Equations: 


m. = ($p $ T + r D Q r D T ) i _ 1 

(3.20) 


/* £ i 



*i = / £ (*i-i> b i . 1 , u i _ 1 , 

Vi 

ti_i)dt + x i _ 1 (3.21) 

“1 

W. = (HMH T + R) i 

(3.22) 

— . 

K v = (MH T W _1 ). 
x i 1 

(3.23) 

_ 

- [I - K x H li S i 

(3.24) 


v i ” y i ' h (*i’ °> u i » 

(3.25) 

n 

*i " + v i 

(3.26) 

n 

i 


$ = 


n 


r D 


1 
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where 


exp {FAT } 
F _1 (4>-I)r 



Bias Filter Equations: 


u i - 4 l-l v i-l * F ^1 fi-l - \ 

K-i 

'(3.27) 

Si - Hi Ui * C t 


(3.28) 

V 1 - U 1 - K Xi s i 


(3.29) 

"i - »i * S i P b 1 . 1 S I 


(3.30) 

T ~-l 

K. = P, S W. x 

b i b i-i 1 1 


(3.31) 

p bi - i 1 - K b s ii ‘'bi.i 


(3.32) 

b i = b i-l * K b. (v i ' S i b i-l 5 


(3.33) 


Note that the primary state estimate at the i t ^ 1 time point 
is denoted by x. and the secondary state (or bias) estimate 

is denoted by b^. The bias filter also utilizes Kx^, and 

Vi which are computed in the bias-free filter. These two 
sets are fully coupled by the composite state update which 
corrects the bias-free state estimate for the effects of 
the bias. 


Composite State and Covariance Update: 
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(3.34) 

(3.35) 

(3.36) 

(3.37) 
47 



For the local smoothing options, it is possible to iterate 
on either the bias-free state estimates and/or the composite 
state filter estimates. 


n 


Bias-Free Local Smoothing: 


x. , | . = x. , + P. , .. [I 

l-l 1 l-l l-l l-l L 




(3.38) _ 


Composite State Local Smoothing: 
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(3.39) 


where 


(3.40) 


The fixed interval smoothing equations are used when the 
filter estimates are inadequate. These equations are neces- 
sary when some states are driven by process noise (i.e., 
gust states) and it is necessary to estimate this process 
noise. A further benefit of fixed interval smoothing is 
that smoothing estimates are computed for all the primary 
states. These equations are shown below. 


, c _ 
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Fixed Interval Smoothing: 
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X i-1 = 

(I - P i H^R: 1 H i ) T [4jA i 
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- hTrT 1 ( z i - H.x.)], A n = 0 

i 

Note 

that x^ 

is the estimate from the bias- free (forward) 



filter. Using the above equations, it is possible to estimate 
the process and measurement noise covariance by: 
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v i v 


T 

i 
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These equations are implemented in the filter/smoother 
algorithm, which is initialized with a priori estimates of 
x , P , Q and R. Figure 3.1 depicts the overall structure of 

the f ilter/smoother , indicating the basic sequence of filter- 
ing, fixed-point smoothing, and fixed interval .smoothing . 

Rotorcraft State Estimation Program . The filter/ 
smoother algorithm of Figure 3.1 is used with different 
rotorcraft models. The reason for some decomposition of 
the rotorcraft model is basically computational; simultane- 
ous quantification of all vehicle dynamic and measurement 
system equations with a Kalman filter is simply beyond the 
memory capability of existing computers (for a given finite 
execution time limitation). The model decomposition is 
based on the principal objectives of an overall rotorcraft 
filtering approach. 

For this effort, the objective is to provide state esti- 
mates and sensor error coefficient estimates for use by sub- 
sequent parameter identification software. The specific 
models (to be included in the filter) which result to achieve 
this objective are as follows: 

(1) Fuselage/gust estimation to estimate the rigid 
body translational and rotational states, errors 
of the associated on-board and off-board sensor 
model, and gust statistics. 

(2) Rotor state estimator to estimate the rotor flap- 
ping and lead-lag states . 

(3) Force and moment static estimator to estimate the 
distribution of forces and moments from load cell 
sensors. (This estimator is denoted as the RSRA 
estimator because its principal application is to 
the Rotor Systems Research Aircraft.) 

Table 3.1 summarizes the resulting filter/smoother options, 
and their applications. 

Figure 3.2 illustrates the basic integration logic for 
the three filter models in order to meet the application 
requirements of Table 3.1. Detailed models for each of these 
options will now be presented. 

Fuselage/gust estimator equations . - The state and mea- 
surement vectors for the fuselage/gust estimator are shown 
in Table 3.2. The user of DEKFIS has the option of choosing 
which states (and, hence, which equations are to be 
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FILTER 

SECTION 


SMOOTHER 

SECTION 


FIGURE 3.1. - STRUCTURE OF DEKFIS. 
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TABLE 3.1. - FILTER/SMOOTHER EQUATION OF MOTION OPTIONS 



FUSELAGE/GUST 

ESTIMATOR 

ROTOR STATE 
ESTIMATOR 

RSRA 

ESTIMATOR 

Single Rotor 
Hel i copter 

/ 

/ 


Multi-rotor 

Helicopter 

/ 

Repeat for 
each rotor 


RSRA-Fixed Wing 
Configuration 

/ 


/ 

RSRA-Compound 

Configuration 

/ 

✓ 

/ 

Whirl Stand 
Data 


/ 



integrated), which measurements and which process noise 
sources are to be used for a particular problem. 


The vector is of size 25x1 and contains the primary 

states for the fuselage/gust estimator. These are: the air- 

craft attitudes — Euler angles <J> , 0, \p; the aircraft iner- 
tial velocities in the aircraft body axis -v , v , v ; the 

x y z 

fuselage angular rates — p, q, r; the linear accelerometer 
indicated acceleration — a X j , ay^ , a Zj (this differs from 


the actual linear acceleration by a first order lag asso- 
ciated with measuring that acceleration) ; the angular accel- 
erometer indicated acceleration — Pj, q^, fj (likewise 

lagged); gust velocities in an inertial north, east, vertical 


and the aircraft position in an 


frame - %- V g E ’ v Sv’ 
inertial north, east, vertical frame — X 


N ! 


T 


The y_ vector is of size 24x1 and contains the measure- 
ments for the fuselage/gust estimator. These are: the air- 

craft attitudes — <b , 0 , ip ; the airspeed in the aircraft 

body axis — v Xm , Vy^, v Zm (or either alternately or redun- 
dantly as V , 8 , a measurements) ; the aircraft angular 
J m m m 
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• FLIGHT DATA 

• PREPROCESSING 



• MODEL STRUCTURE PROGRAM 

• PARAMETER IDENTIFICATION 


FIGURE 3.2. - FLOWCHART AND OPTIONS OF STATE ESTIMATOR 
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rates — p , q , r ; the linear accelerometer measurements — 
r m m 

a XT > ay r , a ZT ; the angular accelerometer measurements — 
i m . /i m > i m 

Pt , Qt , n ; and radar position measurements in an inertial 
r± m ni m i m 

north, east, vertical frame — ^N m > ^E m > Zv m (or alternately 

R "» 8 Rm> “M- 

The x 7 vector is of size 50x1 and contains the biases 
(i.e., b'sj and scale factors (i.e., k's) for each of the 
measurements . 

The w vector is a 9xl process noise vector. This pro- 
cess noise drives the linear accelerations — , w^; the 

angular accelerations — w^, w^, w^; and the gusts — w^, w 7 , . 

The v vector is a 24x1 measurement noise vector each ele- 
ment of wKich corrupts one of the aforementioned measurements. 

The nonlinear fuselage/gust equations are presented below. 

Primary State Equations: 

4> = p + (r cos(J) + q sincji) tan9 

8 = q cos4> - r sin4> 

ip = (r cos<}> + q sin<i>)/cos0 

v x = v y r - v z q - g sine + g sin0(o) + Wj. - w g Jl Zi + w 6 Z yi 
+ (q 2 +r 2 )Z Xl - q p Z yi - r p Z Zl 

Vy = \> z p - v x r + g sintf) cos0 - g sin<}>(o) cos8 (o) + W 2 ~ w^Jl^ 
+ w 4 £ z 2 ' P <1 *x 2 + (P 2 + r2 ) z y 2 ' r * Z z 2 

• v y P + g cost}) COS0 - g COSTCO) COS0(O) 

- w 4 V 3 + Vx 3 - P r £ x 3 - ^ r £ y 3 + 


(P 2+ P 2 H Z3 


v z = v x * 


+ W. 


P = w 4 

q = w 5 

r = \ir 


1 

j 


n 


5 4 


i 



a x : = a l a xi • a l W 1 
a Yl = a 2 a y i ' a 2 w 2 

a z'i = a 3 a z : ’ “3 w 3 

Pj = a 4 pj - a 4 w 4 

«I - o 5 qj - a 5 w 5 

r I = a 6 r I a 6 w 6 

% = “ 7 % " “ 7 “ 7 

\ " a 3 V g E - a 8 w 8 

V §v = a 9 V g V “ a 9 W 9 

Xj^ = cos 9 cos^ + v sin8 sin<}> cosij; - v cos<j> sin^i 

+ v sin0 cosif) costy + v sin<|> sini|j 
Z z 

Yg = cos9 sinijj + sin0 sin<j> sini|> + cos4> cosijj 

+ v sin0 cos4> sin^ - v s in<J> costjj 
z z 

• 

Z,, = -v sin0 + v cos0 sinA + v cos0 cosd> 

V x y Y z Y 

Secondary State Equations: 

*2 = 9. 


Measurement Equations : 

<(> = k, A+b.+n. 

Y m <J> 4> 1 


B = k n 8+b +n 7 

m 0 0 Z 
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k , \ 1 > +b +n 
ip r ip z 


Vx m = 

VvV 

+ 

S + n 4 



V 7m = 

N ( W 

+ 

% + n 5 


— 

v z = 
m 

'‘v^V'V 

+ 

b v z + “ 6 


1 


k p P +b p +n 7 




! 

j 


k q ^ +b q +n 8 




1 

i 

r = 

k T+b +n n 




( 

m 

r r 9 




i 

a y 

X I 

m 

’ NS + 

g 

sin0(o)) + b a 

X I 

* n 10 

j 

a yi 

m 

= k a yi (> yi - 

g 

sin<|)(o)cose (o)) 

+ b a y + n ll 

i 

i 

a " z I 

m 

= ka Zj^ z I 

g 

COS(J) (o) COS0 (o) ) 

* b a z + n 12 
Z I 

"H 

1 

(P^ 

m 

’ k P t Pi + 

b Pl 

+ n 13 


- 1 

I 

l 

(qi) 

m 

= k q t *1 * 

\ 

+ n 14 


1 

| 

( f I>m 

* k r r Z I + 

bz i 

+ n l 5 


t 

X N = 

k x__ X N + b x 

+ 

n 16 


| 

| 

m 

N 

N 




Y E = 
c ra 

k Y E Y E + b Y 

+ 

E 

n 17 


1 

i 

Z V = 
m 

k z v Z V + b z 

+ 

V 

n 18 


-| 
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V m k V ( ( V v x g +q £ z 4 " r l Y ^ 2 + ^y +v 7 g +r £ x 4 "P £ z 4 ) 2 

+ ( V z +V z g + P *y 4 ‘ q ^ 2 Y + b V + n 19 
6 m = k B tan_1 ((V V y g +r £ x 5 “P ^z 5 )/tV v x g +q £ z 5 " r £ y 5 ^) +b B + n 20 

% = k a tan ” 1 ( ( v z +v Zg + P V 6 ‘ q £ x 6 )/(VV q £ V r V) + V n 21 

R m = k R( X N +Y E +Z v) + b R + n 22 

P R m = k 3 R tan'^-Zy/CX^Yg 2 )^) + bg R + 

a R m = k a R tan (" Y £/^^ + b a R + n 2 4 
where , 

v x g = + (V gE +V WE )cos0sin^ - (V gv+ V Wv )sine 

v y g = C v g N + V WN )sin 9 sin<j>cosiJ; - (Vg^+V^) cos<()sini|) 

+ (V +V w ) sin 0 sincj)sin^ + (V +V W )cos<(>cos^ 

°E E °E E 

+ (V gv +V Wv )cos 0 sin<|) 

v_ = (V_ +V, . ) sin 0 cos 4 >cos^ + (V„ +V.. )sin 4 >sinij; 
z g §N W N §N W N 

+ (V. +V W ) sin 0 cos<j)s inili - (V a +V W )sincj)Cos 4 ) 

8 E E g E E 

+ (Vg v +V Wv )cos 9 cos<i) 

The parameters used in these equations are defined in 
the fuselage/gust estimator parameter list shown in Table 3 . 3 . 
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TABLE 3.3. - FUSELAGE/GUST ESTIMATOR- PARAMETER LIST 


ARAMETER 

INDEX 


1-25 

-1 

26-75 

*2 

76 

* X 1 

77 

tyj 

78 


79 

l *2 

80 

hz 

81 

l Z £ 

82 

tx 3 

83 

l y3 

84 

l *3 

85 

* x 4 

86 

l y4 

87 

iz 4 

88 

* x 5 

89 

hs 

90 

l *5 

91 

* x 6 

92 

l y& 

93 

lI 6 

94 

“1 

95 

“2 

96 

“3 

97 

“4 

93 

“5 

99 

“6 

100 

“7 

101 

“8 

102 

“9 

103 

Vw N 

104 

Vw E 

105 

vw v 


PARAMETER DESCRIPTION 


‘s (use uncompressed vector) 
‘s (use uncompressed vector) 


a x accelerometer location, ft (relative to c.g.) 


accelerometer location, ft 


a z accelerometer location, ft 


V probe location, ft 


3 vane location, ft 


a vane location, ft 


accelerometer time constants 


gust power spectral time constants 


steady gust components (wind) 
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Rotor state estimator equations .- The state and measure- 
ment equations for the rotor state estimator are shown in 
Table 3.4. As in the previous estimator, the user has the 
option to flag which states, measurements and process noise 
sources he wants to use. 


The x-l vector is a 14xl vector of primary states used 

in the rotor state estimator. These are: rotor coning angle — 

3 q ; rotor longitudinal flapping — 3^; rotor lateral flapping — 


$ 1S ; similar fixed system coordinates for the rotor lagging 
motion — C Q , £ ^ rotor azimuth angle — and the time 


derivatives of each of these — 3 , 

. o’ 1C’ 


B 1S’ 5 


o’ nc : 


C 1S , and 




R' 


The y vector is 13xi vector of measurements for the rotor 
state estimator. These are: the individual blade flapping 

angles 3i m > i= 1,6; the individual blade lagging angles Ci m> 

i=l,6; and the cosine of rotor azimuth — cos^p m . 

The x^ vector is a 26xl vector of biases and scale fac- 
tors for each of the aforementioned measurements. 

The w vector is a 7xl vector of process noise sources 
each driving one of the seven rotor degrees of freedom 
modeled. 

The v vector is a 13x1 vector of measurement noises each 
corrupting one of the 13 measurements. 

The state and measurement equations for the rotor state 
estimator ' are presented below. 

Primary State Equations: 
x 1 = F ^ + £ w 
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where 


' 1.0 


1.0 


1.0 


F = 


1.0 


1.0 


, r = 


1.0 


1.0 


1.0 


1.0 


1.0 


1.0 


1.0 


1.0 


1.0 


Secondary State Equations: 


*2 = ° 


Measurement Equations: 

6l m = k 0 1 (3 o -6i c cos C^ R + 4> j D - 01 s in x ) ) + b B + n 1 

6i m = k 6 i ( B o" B l c C0S ^ R +( l ) i ) - & l s sin( l p R +c^ i )) + 

B e m = k B 6 ( 6 o" 6 l c C0S (VV -B lg sinC^ R +4> 6 )) + bgg + n & 

= k ?l( ? o“ ? l C C0S ^R +e P -5l s sinC^ep) + b^ + n y 

% = h^o-Zlc C0S( V 6 i } -5i s sinC^Q^) + b 5 . + n i+6 
? 6 m = k S 6 (V ? l C COS( W - ? l s sin( V 9 6 ) ) + b ^6 + n 12 

co sip p = k cosO R ) + b + n - 

m R C0S ^ R 13 

The parameters used in these equations are defined in 
the rotor state estimator parameter list shown in Table 3.5. 



TABLE 3.5. - ROTOR STATE ESTIMATOR PARAMETER LIST 


Parameter 

Description 

$ ^ 


*2 


*3 

Blade flapping measurement 

*4 

phase angles, rad. 

*5 


*6 


9 1 


9 2 


9 3 

Blade lagging measurement 

9 4 

phase angles, rad. 

e 5 


9 6 



RSRA estimator equations .- The state and measurement 
equations for the RSRA state estimator are shown in Table 3.6 

The x-^ vector is a 10x1 vector of primary states for the 

RSRA estimator. These states are: the rotor forces at the 

hub in the aircraft body axis system — X R , Y R , Z R ; the rotor 

moments at the hub — L R , M R , N R ; the inertial accelerations 

of the transmission in the body axis system — a t , a t , a t ; 

x y z 

and the sum of the engine and tail rotor drive shaft torques 

V 


The y vector is a 10x1 vector of measurements consisting 
of: the transmission load cell reactive forces — A, B, C, D, 

E, and F; transmission accelerations — a t , a t , a-t > anc * 

x m ^m z m 

total drive shaft torques — Q . 

The x ^2 vector is a 20x1 vector of biases 
tors for each of the 10 measurements. 


and scale fac- 








The u vector is a 6xl control vector of body angular 
accelerations — p, q, r; and body angular rates — p, q, r 
that are treated as deterministic inputs. These six quan- 
tities should be obtained from a prior fuselage/gust esti- 
mator run. 

The w vector is a 10x1 vector of process noise sources 
that drives each of the states. 

The v vector is a 10x1 measurement noise vector. 

The state and measurement equations for the RSRA state 
estimator are presented below. These equations are based on 
Reference (11). 

Primary State Equations: 

Xi = w 

Secondary State Equations: 

k 2 = 2 . 


Measurement Equations: 

A ■ h U X R +h 12 Y R +h 13 Z R* h 14 L R* h lS M R +h 16‘\' h ll m t at x 
-h 12 m t a tjr -h 13 m t a tz -h l4 It xx P-h 14 (I t ^-I t ^)r p 

* h 14 £ 6 m t at y" h 14 C *t' h 15 It rr q ' h 15 l ' It xx’ It zz' lq r 
+h 15 f 5 m t a tz- h lS f 6 m t a tx‘ h 16 I tzz f 
■ h 16 f 5 m t a ty‘ h 16 tlt yy' lt xx )pq * "l 
B " h 21 X R +h 22 Y R* h 23 Z R +h 24 L R* ll 25 M R* ,I 26 N R' h 21 m t a t x 

- h 2 2 '"t a t y - h 23 m t a t 2 - h 24 I txx P ' h 24 CI tzz‘ It yy )r p 
* h 24 f 6" , t a V h 24Qt- h 25 I t rr q - h 25 (I txx' I tzz )q r 


64 



+h 2S f 5 ra t a t z - h 2sV't a V h 26 I t 22 f 

- h 26 f 5 m t a V h 26 (I t yy - I t xx )P1 + n 2 


C * T7^T X R ’ (T“fTJT Z R ‘ (T^] M R 


(Tpip”t a 'x 


C! t -I t ) 

L xx L zz 


f 5 m t 


l yy • 

+ (f 1 + f 2 ) m t a t z + Cf 1 +f 2 ) q + Cf^f 2 ) qr " (f 1 + £ 2 ) a t z 


f 6 m t 


+ (T^ a t x + n 3 


D = 


<Ty^l Y R " (£-^ + £ 3 J + (f 1 + f 3 3 m t at y + (±^+±” 3 ) 


d t - x t ) 

L yy L xx 

(£i + £ 3 ) 


pq 


f 5 m t 


(•fi + £ 3 ) 


a t„ + 


f l +£ 3" £ 2 
E = i r £ Y 

£i +f 7 R 


■N 


f l +f 3' £ 2 


l + £ 3 r ( £ i + f 3 ) R £ l + £ 3 1 Ly C £ i "*■ £ 3 ]) 


a t. 


'ZZ 


d t - : t ) 

yy xx 

(f 1+ £ 3 ) 


pq 


f 5 m t 


a t + n c 

( £ i +£ 3 ) / 5 


F = X D -m. at- + n 

K P X 


a t x = a t x + n 7 
x m X / 


a t = a t + n Q 

y 8 
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a t 


= a t + 

z m 2 


n. 


% - Q t + n io 

where , 


l ll 

h 21 

2 (fx+fa) 



- f 3 , <V f 3> f 2 

12 

" h 22 

yt yt( f i +f 3) 


= h = 

1 f 2 

13 

23 

2 Cfi+f 2 ) 

14 


i 

-h 24 



— V\ — 

i 

15 

n 25 

2(fl+f 2 ) 


= -h = 

CV f 3 ) 

16 

n 26 

y t (£i + f 3 ) 


The parameters used in these equations are defined in 
Table 3.7. Figure 3.3 shows the RSRA transmission orientation 
on the airframe and the location of the six load cells (A-F) . 
Figure 3.4 shows the physical description of the parameters. 
The most basic parameters are x t , y t , z^, z R and i where: 

x t ’ ^t ’ z t are rotor counting geometry dimensions in the 

transmission principle axis system; z^ is the distance from 

the rotor hub to the transmission center of gravity along the 
shaft; and i is the transmission incidence with respect to 

the longitudinal body axis. 
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Parameter 


Description 


y T 

See Text 

f i 

f , = T cos *t + Cz R + z t> 5ln 't 

f 2 

f 2 = T cos i t “ tz R + z t } S1 ' n ] t 

f 3 

x t 

^3 = ^ Z R + z 0 cos H + ~T s ^ n H 

f 4 

x t 

f 4 = Cz R + z t } cos i t “ T Sin 1- t 

f 5 

f 5 = Z R sin 't 

f 6 

f 6 = z R cos i t 

M t 

Transmission mass, slugs 

*txx 

l tyy 

*tzz 

^Transmission principle moments of inertia, 
(slug - ft2. 







FIGURE 3.3.- RSRA HUB FORCE/MOMENT MEASURE- 
MENT SYSTEM CONFIGURATION 
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FIGURE 3.4.- RSRA ROTOR HUB FORCE/MOMENT MEASURE- 
MENT SYSTEM PARAMETER DEFINITIONS 
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CHAPTER IV 


MODEL STRUCTURE ESTIMATION 


Requirements 

System identification methodology has been used success- 
fully in many applications. This wide scope of application 
has focussed attention to one particular phase of the system 
identification methodology- - the model structure determination 
phase--as an essential step. This step consists of proces- 
sing the input/output data to determine the significant 
linear and nonlinear equations and associated parameters 
which are necessary to represent an observed system response. 
The significance of model structure determination in the 
overall system identification procedure is illustrated in 
Figure 4.1. A number of techniques must be used in the data 
processing stage, particularly for nonlinear operation re- 
gimes, to obtain the maximum information from the data for 
any particular rotorcraft. In the first stage, unmeasured 
or failed channels of data are reconstructed based on avail- 
able measurements. This also gives preliminary force and 
moment coefficient time histories of interest. In the model 



FIGURE 4.1.- ROLE OF MODEL STRUCTURE DETERMINATION IN 
ROTORCRAFT SYSTEM IDENTIFICATION PROCEDURE. 







structure determination (MSD) stage, the dependent forces 
and moments are related to the independent variables (e.g., 
angle-of-attack, Mach no., etc.) to provide the model with 
the most useful predictive capability. The maximum likeli- 
hood method is used in the third stage to refine the para- 
meters of the model structure selected in the previous stage, 
as well as to validate the model structure. 

It must be noted that this procedure will produce a 
mathematical model from the data which represents the par- 
ticular data characteristics. No one set of data can com- 
pletely define a unique model for all possible conditions. 
Additional dynamic tests (with specified input designs) , 
scale model tests, or analytical correlation are always 
required to arrive at such a unique model. These additional 
requirements are discussed in detail in this chapter. 

This step is particularly important in rotorcraft model 
identification because: (1) general rotorcraft models tend 

to be of high order with many dynamic effects, (2) rotor 
dynamics may or may not be important in any maneuver, (3) 
longitudinal and lateral motions may be coupled, and (4) 
there are significant nonlinearities in aerodynamic behavior 
particularly near hover and in transition. 

Clearly the model structure determination is an extreme- 
ly important step in rotorcraft identification work. Yet, 
the authors think that this is the first attempt at MSD for 
rotorcraft. The basic rotorcraft characteristics add new 
complexities to the problem. This chapter will also discuss 
how some techniques previously used for fixed-wing airplanes 
are specialized to rotorcraft applications. Numerical ap- 
proaches are also described because the number of models to 
be tested may be very large. 

The second section discusses the general model struc- 
ture estimation method. Specializations required for rotor- 
craft applications are given in the third section. Finally, 
the fourth section describes the basis of a specific algo- 
rithm on which the computer program developed in this effort 
is based. 


Model Structure Estimation Approaches 

Though the fundamental problem of model structure esti- 
mation has been recognized since the origin of physical 
science, the utilization of statistical methods to systemat- 
ically determine significant characteristics in nonlinear 
dynamic systems was first presented in 1974 [12,13]. Such 
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methods have been applied to static systems [14] (e.g., 
in econometrics and biometrics) and to linear systems 
[15]. From a more fundamental approach, however, the 
work of Akaike [16] and Kullback [17] provided the sig- 
nificant theoretical formulation concepts for both linear 
and nonlinear systems. The theoretical basis of the methods 
discussed in this chapter is found in the cited work of 
Thiel, Akaike and Kullback. The application to aerospace 
systems can be traced to the work of Gerlach [18] in the 
Netherlands who demonstrated the use of a linear technique 
on actual test data. 

The model structural determination process for rotor- 
craft has been isolated into three significant issues: 

(1) Specification of classes of a priori models. 

(2) Criteria for hypothesis testing of data against 
models. 

(3) Numerical procedures. 

Specification of classes of a priori models .- The 
first major step in the model structure estimation process 
is the selection of mathematical forms to be correlated and 
tested against data. The specification of these models must 
be sufficiently broad to include the most probable relation- 
ship without attempting to consider all models. In order 
to form a tractable a priori model base, the following gen- 
eral considerations have been found to be relevant: 

(1) It is useful to recognize two aspects of mathemati- 
cal models of physical systems. Basic models are 
derived from fundamental physical laws. Useful 
engineering models are, to some extent, empirical, 
wherein many of the system interrelationships are 
defined by tests and appropriate interpretation of 
the results. This obviously means that the role 

of physical engineering analysis in the model 
structure estimation process is essential. In 
particular, it is important to keep in mind that 
rotorcraft may have faster and slower modes. 

(2) The ultimate validation criteria, which include 
ability to explain existing data, consistency of 
empirical results with phenomenological considera- 
tions and ability to predict data, dictate model 
form. 
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(3) The objective for which the model is to be used 
will indicate the levels of complexity required 
and the regime over which the model should be 
valid. This means that there may not be a unique 
model to describe a particular vehicle. In fact, 
model selection depends mainly on the use to which 
the model is put. 

Selection of Functional Forms.- In order to optimize 
the initial model for computational speed and numerical 
accuracy, it must include elements of known phenomenological 
behavior supplemented by various classes of functional 
forms, depending on the particular problem. In the absence 
of strong phenomenological evidence to the contrary, poly- 
nomials are well known as simple and versatile choices for 
mathematical models. Indeed, polynomials serve exceedingly 
well in modeling lift, drag, and moment coefficients for 
helicopters. If, however, it is known that the process is 
susceptible to a growth phenomenon, such as vortex buildup 
from aircraft control deflections, exponential or logarithmic 
functions are desirable. The use of polynomials should 
therefore be undertaken with a reasonably clear understanding 
of the limitations. This fact was emphasized in the second 
chapter. Factors which affect the particular polynomial form 
including the following: 

Cl) A priori phenomenological information about the 
process . 

(2) Range of variation of the dependent variables (and 
possibility of differences in physical phenomena 
over that range) . 

(3) Specification of whether the polynomial is to be 
differentiated (requiring valid slope representa- 
tion) or integrated (requiring high reliability in 
the initial value) . 

(4) Computational resources available to use the model 
(in terms of speed and memory). 

To meet these requirements, two basic polynomial formu- 
lations have been used, as shown in Table 4.1: 

(1) A regular polynomial is used to represent a simple 
continuous phenomenon. This polynomial is usually 
linear in the unknown parameters (though many 

cases arise in which the polynomial is nonlinear 
in the unknown parameters) . For improved numeri- 
cal conditioning, orthogonal polynomials are 
desired. The classical orthogonal polynomials of 
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TYPE OF 
PHENOMENON 


BASIC FORM 


Homogeneous 

POLYNOMIAL 

m . 

y = l C.x 
i=0 1 


SPLINE 

m . n m 

y = E C.x + E E C. .(x-x.) 1 

1-0 1 j=l i=v+l 1J J 

l 0 x < x. 

where (x-x . )' - < . J 

J * (x-x ■)' X >_ X. 

J J 

Heterogeneous 














Legendre, Laguerre, and Hermite may be generated 
by three term recurrence relations which are easily 
programmed. The Tschebycheff polynomial is known 
to demonstrate the properties of both the Fourier 
series and the orthogonal polynomials [19]. 

(2) A spline of order n and continuity v is used 
to represent a simple heterogeneous phenomenon. 

The spline function consists of piecewise poly- 
nomials wherein the derivatives are continuous. 
Hence, the spline preserves the continuity of 
lower order derivatives across function discontin- 
uities. There. are n+1 regions defined by 

x < x. , x. < x < x- , . . . ,x < x. At each of the 
11— i n — 

transition points, x^, which are called knots, 

the first derivatives of the function y are 
continuous. The independent variable y is a 
linear function of parameters, C^, but is a non- 
linear function of knot locations. B-splines are 
used for improved numerical conditioning (see 
Table 4.1) . 

The representations for a single independent variable 
can be generalized to many independent variables in several 
ways. Straightforward generalizations of the polynomial and 
spline forms for two variables are given in Table 4.2. Such 
an approach is useful for cases where the number of indepen- 
dent variables and/or terms in approximating polynomials is 
small. For several independent variables, a more organized 
procedure has been found necessary. The method which has 
been implemented to achieve this is an extension of 
Ivakhnenko's group concept [20]. Table 4.2 shows the typi- 
cal equations which are used to evaluate large levels of 
subsets of model terms. For example, the equation of Table 
4.2 (for regular polynomials) is written in terms of the 
second variable for all the unknown variables in the first 
equation. 

The models specified using the above procedure will 
in general be too complex. The next section discusses sta- 
tistical criteria which are used to simplify the above rep- 
resentations such that useful models are obtained from a 
limited set of data. 

Quantitative criteria for comparison of competitive 
models - against test data .- The representations of the pre- 
vious sections are quite flexible and allow for a very large 
number of model structures. In certain applciations , several 
competing representations based on different independent 
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variables may be hypothesized leading to a further increase 
in the number of plausible model structures. The use of 
measured data isolates the most likely model and indicates 
model adequacy in explaining the observed behavior. This 
section discusses several classes of criteria used for this 
purpose. The tradeoffs to be considered in the particular 
selection of the criteria are: 

Cl) The distribution function of the noise. Certain 
criteria are applicable when the noise is white 
Guassian while others are more general; 

(2) A priori knowledge about noise distribution func- 
tion; and 

(3) Number of models to be compared. 

Let there be N sets of measurements (or reconstructed 
values) represented by y^ and x^l), x.^2), ...» x i (p), 

i = l,2,...,N. Quantitative criteria used for model substan- 
tiation based on this data may be divided into four broad 
categories, shown in Table 4.3. Note that all these criteria 
can be used with both the equation error and the dynamic 
model formulations. 

Fit error statistics. Fit error is a measure of the 


difference between the measured response 

and 

its estimate 

based on the 

model. Suppose 

two models 

M 1 

(with 

m^ para- 

meters, 0^) 

and M 2 (with 

m 2 parameters, 

0 2 ) 

are to 

be compared. 

Let y^ (M x , 0 ) 

and y i (M 2 , 

§ 2 ) 

be the 

esti- 

mated values 

A 

of y. based on 

models M^ 

and 

m 2 . 

respec- 


A A * 

tively (0-^ and 0 2 are the corresponding parameter esti- 


mates). Then the fit error leads to the following criterion 


* 1 Select Model 

(4.1) 

> 1 Select Model M 2 


Improved results are obtained if the fit error is corrected 
for the number of unknown parameters in the model. The 
adjusted fit error for model M^ is 


ft . E [yi-yiWi.Spi 2 

1=1 
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TABLE 4.3.- COMPARISON OF CRITERIA FOR VALIDATION OF 
MODELS AGAINST TEST DATA 


CLASS 

EQUATION 

COMMENTS 

FIT ERROR 

• Error Covariar.ca 


• Since fit error always increases with 
number of parameters, subjective termination 
criteria required. 

• Whiteness Test 

N ^I'^i^i+l'^l+l^etc. 


• Fit Error Corrected for 
Degrees of Freedom 



LIKELIHOOD APPROACH 
• Like] ihood Ratio 

p(y|m 1 .5 1 j/p(y m 2 .§ 2 ) 

• Equations given for two models—see text for 
generalizations 

• Requires knowledge of probability distribu- 
tions, a priori 

• Log Likelihood Ratio 
Corrected for Degrees 
of Freedom 

log{p(Y|M 1 ,3 1 )/p(Y|M 2 ,3 2 } 
-2irij + 2m 2 

• Works with any noise distribution 

PREDICTION ERROR 



• Direct Determination 
over an Independent" 
Data Set 

1 Np 2 

• Excellent when model is used for prediction 

• Estimate Over Same 
Region as Data 

N+m, , N , 

• Automatically incorporates degrees of freedom 

• Estimate Over a New 
Region 

N n 

1 P ? 

I r z [af +af(M-)] 
P i-1 1 1 


F-RATIO 



• Equation F-Ratio 

R /m 

( 1- R 2 ) / ( N-m ) 

• Most used test in econometrics and biometrics 

• Parameter F-Ratio 

(R 2 -R 2 )/m 2 

(l-R 2 )/(N-m r m 2 ) 

• Easy to implement 

• Parzen' s Test 

See Ref. 13 

• Excellent for a first cut 
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(4.2) 


1 

(N - m 1 ) 


N . 2 

l [y i -y i (M 1 ,0 1 )] z 


The fit error should be rarely used directly for comparison 
of models. It is, however, very useful in establishing the 
validity of the model selected by other approaches. 

Likelihood ratio statistics. The likelihood approach 
has been used extensively for both parameter estimation and 
comparison of competitive models based on test data. The 
central concept is the likelihood function, which defines 
the probability that the measured data were generated by any 
specific model or any set of parameter values. The likeli- 
hood function for model is p(Y|M^,0^) where Y is 

the set of measurements y^y^,...^^ and 0^ is the maxi- 
mum likelihood estimate of 6^ assuming model holds. 

The model selection criterion is the ratio of the likelihood 
functions for models and M 2 * The likelihood ratio 

must, however, be corrected for the degrees of freedom. 
Otherwise the results will always favor more complex models. 
When there are p competitive models to be compared, the 
following procedure may be used. 

Step 1: Compute J^, i=l,2,...,p 

J i = In {p CY | M i , 6 ± ) } - 2IIU 


Step 2: Let J^. be such that 


J. > J. , i 5* j 


Then is the "most likely" model. 

For two candidate models the method simplifies to 


( p(Y |M, ,0. ) ) / >0, Select 

In { - 1 _ 2m-,+2m 7 > 1 

(p(Y|M ,0 2 )J 1 | <0, Select M 2 


(4.3) 


Important generalizations of the likelihood ratio result 
when one model is the subset of the other model, or in gener- 
al, all models considered are subsets of the same maximum 
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model. Let the maximal model be characterized by parameters 
0. When certain parameters are zero, a lower order subset 
is obtained. (Note that this may be used to model a wide 
variety of situations.) The log likelihood function may 
then be expanded in^a multidimensional Taylor series about 
our best estimate 0 of 0. 


In p (Y | 0 ) = Jin p (Y j 0) + - 3 - P 1 ^ (e - 0) (4.4) 

30 

+ \ (0 - 0) T M(0 - 0) + ... 


where M is the information matrix. Since 0 is the best 
estimate of 0, the first gradient of the likelihood has 

mean zero and covariance M~^. To test the hypothesis that 
a lower order model is valid, the estimated 0 equal to 
zero will decrease the likelihood function significantly. 

A model structure determination principle based on this con- 
cept was detailed by Gupta [21]. 

The likelihood method is optimal under a variety of cir- 
cumstances. Its rigorous applicability is limited by the 
theoretical requirement to know the probability density of 
the measurements for each model. 


Prediction error statistic. - The capability of a model 
to predict system responses for a class of inputs is a 
desirable quality. Therefore, prediction error of models 
may be used as the quantitative criteria to select the model 
which best substantiates the measured data. There are two 
methods to determine the prediction error of a model. In 
the first method the measured data is divided into two parts. 
The first part of the data is used to estimate unknown para- 
meters in each model. The estimated models are then used 
to predict the response for the second data set. The model 
which predicts the response most accurately is the desired 
model. The second method is indirect where the prediction 
error is estimated statistically. The prediction error is 
either based on the same input as the one used to estimate 
the models or covers a different region (the choice depends 
on the ultimate application of the model). A good estimate 
of prediction error for a model with m^ parameters 

over the data region may be shown to be [16] 


N+m^ 


£ [y.: - yiCM.,6,)] 
i=l 1 1 1 1 


(4.5) 
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The prediction over a data region y ^ 
computed as follows. Let 
value for y^ with mean square error 


M. 


1 * 


If 


a . 

1 


is the variance of noise 


tion error is 


, i = 1 , 2 , . . . ,N is 

p 

be the estimated 

/v 2 

for model 

in y .., the predic- 


N P 

z + ,)] • (4.6) 

p i=l 1 11 

Prediction error criteria are most useful when the estimated 
model is to be used for simulation. They automatically in- 
corporate the degree of freedom correction for number of un- 
known parameters in the model. 

F-ratio statistic.- The F-ratio is perhaps the most 
widely used statistic for model hypothesis testing 
in econometrics and biometrics. The test is based on 
the assumption of normally distributed random disturbances 
and requires a priori specification of acceptance-rejection 
boundary. 

These assumptions are restrictive, and as noted in Ref. 
12, F-tests should not be used as the only criterion for the 
adequacy of a model. On the other hand, they do have com- 
pensating attributes which include the following: 

(1) Standard algorithms have been optimized for com- 
puter implementation of F-ratio statistical hypo- 
thesis testing. They allow the consideration of 
an extremely large number of models. 

(2) A relative maximum of the F-ratio with the number 
of parameters is often found in practice. This 
maximum is cited in Ref. 12 as an experimental 
result, which to the authors' knowledge, has not 
yet been investigated by other researchers. 

(3) For many practical cases, it gives the same result 
as more sophisticated approaches. 

The desirable performance of the F-ratio, in spite of the 
assumptions, is presumably due to the robustness of the 
statistic and the particular self-check features used in the 
implementation. In particular, application of the F-test 
criterion on properly prefiltered data has been found to 
improve the utility of the statistics. Three approaches 
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have been used based on: (1) equation F-ratio, (2) parameter 

F- ratio, and (3) Parzen's test. The equation F-ratio tests 
the validity of the entire model. It is given by 


R 2 /m 

(1 - R 2 )/(N-m) 


(4.7) 


where R is the equation multiple correlation coefficient, 

N is the number of data points, and m is the number of 
parameters in the equation. 

The parameter F-ratio statistic is applied as follows. 
Suppose that the complete set of parameters 0 is divided 
into two subsets of 8^ and 02 of size m^ and n^. An 

F-test can be used to test the hypothesis that ©2 is equal 

to any specific value (usually zero) while parameters 0^ 

are chosen to minimize fit error. The F-ratio of parameters 

0 2 is 

(R 2 - R?)/m 7 

F(0 ? ) = ^ ± * (4.8) 

(1 - R^)/(N - m 1 - m 2 ) 

where 

R, is multiple correlation coefficient with para- 
meters 0^ 

R 2 is multiple correlation coefficient with para- 
meters 0^ and ©2 

If the F-ratio is small compared to a threshold, parameters 
©2 may be set to zero. F-ratio for parameters in the equa- 
tion is computed in a similar manner. In practice, the F- 
test is performed on single parameters rather than on sets 
of parameters. 

Parzen [22] devised a unique method based on the F- 
statistic to develop models for physical systems. The 
method assumes that the "true" model is very complex and 
has several degrees of freedom. Since the measured data 
contains noise and does not encompass the entire operation 
regime, the estimated model must be a simplified version of 
the true model. Each simplified model is then compared with 
a high order model using an F-statistic. 
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Summary.- Quantitative criteria presented in this sec- 
tion look at model structure optimality from different view- 
points. Extensive experience has indicated that for well- 
behaved systems, all the tests (except fit error) give simi- 
lar model structures if the noise has a normal distribution 
and the specified a priori minimal model includes all the 
effects observed in the data. If the distribution of the 
noise is significantly different from normal or if the noise 
is not white, the likelihood ratio test gives the best 
results . 

Numerical procedures .- Implementation of methods to 
develop models from test -data must consider the following 
factors : 

(1) The number of hypothesized models may be very 
large (number in billions or more is common). 

(2) All models contain some level of modeling error 
(i.e., no model of a physical system is perfect). 
The models are "good" or "poor," not "right" or 
"wrong . " 

(3) Many physical systems are dynamic. The modeling 
of system dynamics may improve the model structure 
estimates . 

(4) Distribution functions are usually not known and 
must be approximated. 

(5) It should be possible to incorporate the analyst's 
opinion. 

(6) Computation time should be reasonable . 


These factors indicate that proper implementation is a key 
part of the successful model structure determination process. 

Three formulations have been successfully used: 

(1) Equation error. 

(2) Kalman filter (or Extended Kalman filter) . 

(3) Maximum likelihood (or its special case, output 
error) . 

Equation error formulation.- In the equation error for- 
mulation, measurements of dependent and independent variables 
are related by 
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(4.9) 


y ± = £(X i ,0) + e ± , i = 1,2,. . . ,N 

X i is the set of p independent variables at the ith 
point and 0 is the kxl vector of unknown parameters. 

and X^ are obtained by reconstruction from available 
measurements. For example y^ may represent angular accel- 
erations obtained by differentiating measurements of angular 
rates. Such reconstruction is often necessary in equation 
error formulations. The results shown here are for uncor- 
related and Gaussian noise, though extension to other cases 
is straightforward. The above equation is linearized about 
the nominal value 0 of the parameter vector 


9 f (X • 0 ) 

<■* - V * H 

0 

i = 1,2, ... ,N 

which may be written as 

AY = AA9 + e (4.11) 

AY is an Nxl vector and A is an Nxk matrix. Note 
that Eq. (4.1) will usually be an overdetermined system 

(N > K) . A0 is chosen to minimize | | AY - AA 0 | | ^ . 

Nominal parameter values 9 0 are obtained either from 

a priori estimates or from a previous iteration. Equation 
(4.11) must be solved a number of times until convergence 
occurs. Most nonlinear equations of engineering signifi- 
cance can be solved by iteratively forming and solving a 
series of linear equations (4.11). 

A straightforward way to evaluate all of the possible 
models is to solve the complete least square problem for 
each of the parameter subsets. This is unpractical even 
for systems with 25 unknown parameters (34 million models). 

One method to make the procedure feasible is to use the 
stepwise regression procedure [12]. By considering one 
parameter at a time, only a small fraction of all the sub- 
sets is tested. This allows the analysis of models having 
up to 400 candidate parameters. One disadvantage of the 
method is that it usually finds only one subset of each size, 
unlike the complete search methods. 
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Kalman filter.- The Kalman filter and extended Kalman 
filter have been applied to many problems in state estima- 
tion, parameter identification and fault detection. Like 
the least squares method, the parameter estimates resulting 
from a Kalman filter are biased. To use this approach for 
model structure estimation, one set of dynamic equations is 
written for each model proposed, i.e., for the ith model 

x. = f(M i ,x i ) + w(t) , 0 < t < T (4.12) 

y(k) = h(M^,x^,k) + v(k) , k = 1,2, ...,N (4.13) 

The states x^, of course, include the unknown parameters 

in each model. An extended Kalman filter is developed for 
these equations. One or more of the following quantities is 
then used to compare the models: 

(1) Innovations : The innovation for model i is 

A 

y(k) - y (k, k- 1 ,M^) . Its bias and covariance is 

a good measure of one-step-ahead prediction error. 

It could also be used to compute certain likelihood 
ratios by making suitable assumptions. The prob- 
lem with innovation is that while the parameters 
are being adjusted in the initial portion of data, 
innovations are large and more complex models may 
be unnecessarily penalized. 

(2) Fit Error: The fit error for model i is 

y (k) - y(k N,M^) and can be used like innovations 
to compare models. Residuals may also be used. 

(3) Parameter Estimates and Covariances : This is a 

poor basis of a test because predicted covariances 
are often grossly in error when determined using a 
Kalman filter. 

Problems with using a Kalman filter for model structure 
determination are: (1) choice of measurement and process 

noise covariance have a strong influence on model; in general, 
increasing process noise covariance will result in less com- 
plex models, (2) filter divergence because of poor starting 
values may invalidate an otherwise good model, and (3) since 
one extended Kalman filter is required for each model (34 
million Kalman filters for a problem with 25 parameters, 
each of which could be zero), the method often requires un- 
acceptable computation time in practical systems. 
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Maximum likelihood method.- Several of the problems 
associated with extended Kalman filter may be solved by 
using the maximum likelihood method. This method has been 
described extensively [12,13] in the literature. It gives 
accurate estimates of parameters as well as associated error 
covariances. Though all of the quantities used for model 
comparison with Kalman filter may be applied, the liklihood 
ratio and prediction error tests are most appropriate. In 
certain cases, log likelihood ratio expansions of the kind 
shown in the likelihood statistic test section can simplify 
the implementation significantly. 

The maximum likelihood procedure is very general and 
includes equation error and Kalman filter procedure as spe- 
cial cases. A major problem with a direct application of the 
procedure for model structure estimation is the excessive 
computation time requirement. This technique is best applied 
when the equation error has substantially reduced the number 
of plausible models. 

Summary . - A comprehensive advanced model structure 
determination for dynamic systems is an essential aspect of 
system identification for rotorcraft. The critical elements 
of such a procedure are: 

(1) Specification of functional forms for useful 
models . 

(2) Criteria for selecting an adequate model against 
competing models. 

(3) Efficient numerical algorithms for integrating the 
conflicting requirements of adequate functional 
forms and accurate criteria. 

The next section will show how these choices are made 
for rotorcraft model structure determination problem. 


Rotorcraft Model Structure Estimation Method 

The previous section discussed a fairly comprehensive 
model structure determination method for a wide variety of 
rotorcraft. Conceptually, any combination of functional 
form, model evaluation criteria and numerical procedure may 
be used in developing mathematical models of rotorcraft from 
real time data. Practically, the choice of an acceptable 
combination is dictated by the following considerations: 
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(1) Rotorcraft models are usually of high order with 
complex interactions, both mechanical and aero- 
dynamic between rotor(s) and fuselage. 

(2) There are many different classes of rotorcraft and 
very little experience exists in processing rotor- 
craft data. 

(3) There is poor wind tunnel data base for many 
helicopter types and analytical simulations (like 
C-81) though useful for simulations are not as 
good for determination of model form. 

(4) The characteristics of typical rotorcraft change 
significantly from hover to transition and cruise. 
In addition, functional relationships between forms 
of aerodynamic coefficients at various speeds is 
not thoroughly understood. 

(5) Because of model complexities, certain important 
interactions are likely to be unidentifiable for 
any maneuver. 

(6) In the face of the above problems, the computation 
time should be reasonable even for complex and long 
maneuvers in which several control inputs are 
applied sequentially or simultaneously. 

These considerations suggest the following selections: 

(1) A model based on dynamic equations of motions, in 
which all unknown functional forms (aerodynamic 
coefficients, interaction coefficients, unknown 
torques, moments and forces) are written as multi- 
dimensional polynomials in corresponding indepen- 
dent variables. 

(2) F-ratio, equation F-ratio or prediction error 
criteria are used in model selection. 

(3) The equation error method is used in the model. 

This is partly justified because the MSD step is 
followed by a parameter identification step based 
on the maximum likelihood technique. 

We now describe an algorithm based on the above selec- 
tions in detail. 



Description of a Model 
Structure Estimation Algorithm 

The use of the state estimation approach as the first 
step on the data, and the selections of the polynomial func- 
tion forms and equation error approach lead to a model struc- 
ture problem in the following general form 

y = X0 + e (4.14) 

where y is an (mxl) vector of observations, X is an mxp 
matrix of known coefficients (ra > p) and 0 is a pxl vector 
of parameters. e is an (nxl) vector of disturbance or error 
random variables with = 0 (zero mean) and variance 

T 2 

^(ee ) = a I (components of £ are uncorrelated with the 
same variance a 2 ). 

Basic principles .- Minimization of J with respect to 0 
yields the well-known least squares estimator 


0 = (X T X)' 1 X T y 


T 

where (X X) is nonsingular. Substituting Eq. (4.14) into 
Eq. (4.15), it is easily shown that 


<?(0) = 0 . (4.16) 

<?[(0-0) (0-0) T ] = a 2 (X T X) _1 (4.17) 


where X is fixed in repeated samples and assuming 
<F(X T £ ) = 0. 

? 

In actual experimentation, a is not known and must 
be estimated. This is due to the fact that e = y -X0, the 
error, is a nonobservable, stochastic variable. However, 

estimates of the variance of e = y - X0 (Fig. 4.2), s , may 

2 

be determined which estimate 0 . The sum of squares of 
the residuals is found by solving Eqs. (4.11) and (4.14) for 

e c, i . e • , 


aTa ^ T 1 T* 

£ 1 e = (y - X0) 1 (y - X0) = e 

= y T ~tfy 
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where ^t( = (I - XOX T X) _1 X T ) . Note that 2 =^U . Since 

is thus idempotent and tr«i/ = m-p,* taking the expec- 
tation of Eq. (4.18) yields 


5 2 = 5 ^ Cy - X0) T (y - X0) 


aI* /V 

e e 
m-p 


(4.19) 


2 

which estimates a . Equation (4.19) in Eq. (4.17) gives 
the sample estimate of the covariance of § 


(0 - 0) (0 - 0) T = s 2 (X T X)" 1 


(4.20) 


Note that Eq. (4.15) produces unbiased estimates of 0 ohly 
if the model is correct. 


Decomposition of variance .- Substituting the estimate 
into J, it is found that 


J 


aTa T aT a 

ee = yy-yy 


(4.21) 


where 

y = E (y | X) = X0 . (4>2 

Equation (4.21), a consequence of orthogonality theorem of 
least squares, decomposes the sum of squares of the data 
T 

variation (y y) into the contribution from the sum of 
squares of the regression equation (y y) and sum of squares 

/N'T 1 a 

of the residuals (e e), i.e., 


T aT a t 1 a 
y y = y y + e e 


(4.23) 


The basic idea of subset regression is to compute the 

a'T’/n m 

reduction in error sum of squares (e e) relative to the 

rp 

increase in regression sum of squares (y y) which is 


*The rank of I - = X(X^X) ^ is p, and the rank of I 

is m. If the observations are not centered, then p is 
replaced by p+1. 
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caused by adding (or subtracting) new variables 0^^ in the 

regression. Ratios based on these incremental sum of squares 
may then be used to determine the significance of these 
added parameters. These ratios may be formulated as test 
statistics for which, under the assumption of normally dis- 
tributed errors, e, standard significance testing may be 
performed. That the errors are, in fact, normally distri- 
buted, is justified by the central limit theorem on the nor- 
mality of the distribution of a large number of random dis- 
turbances . 

A 

Significance of the regression equation .— The solution 0, 
Eq. C 47 1 5 ) , to the minimization of error sum of squares is 
unique when X has full column rank. Any coefficient vector 
which differs from the least squares vector 0 must lead to 
a larger sum of squares. That is, for any other coefficient 
vector, 0, 


(y - X0) = (y - X0) - X (0 - 0 ) (4.24) 

has an error sum of squares 

(y-X0) T (y-X0) = y T y + ( 0 - 0 ) T X T X( 0 - 0 ) . (4.25) 

A 

The quantity (0-8) is the "sampling error" of the regres- 
sion. 


Assume that it is desired to test the hypothesis that 
0 is some specific numerical value, say 0 Q (e.g., 0 Q = 0 

if there is not dependence of the data on 0). Then, if 

0 is the true value of 0, 
o 

0-0 Q = 0-8 = (X T X)' T X T (X0+e) -0 = (X T X)' 1 X T e. (4.26) 

Using Eqs . (4.18) and (4.16), Eq. (4.23) then becomes 

(y-X0 o ) T (y-X0 o ) = e T e + e T (I . (4.27) 

In order to fulfill the requirements for using statis- 
tical tests based on the sum of squares of Eq. (4.27), it is 
necessary to recall the fact that an indempotent quadratic 
form in independent standardized normal variates is a chi- 
squared variate with degrees of freedom given by the rank of 
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the quadratic form* (Sec. 15.11, Ref. 23). Furthermore, if 

2 2 

the two normalized sum of squares S./a and S./a are 

2 2 1 J 

X (v.) and x (v ■ ) distributed, and S^ and S^ are 

independent, then 


F 


V v i 

ST7^ 


(4.28) 


is a Fisher F-ratio distributed with v. and v. degrees 
of freedom. 1 J 

As discussed after Eq . (4.18), the rank of is m-p 
and that of I -«✓, U is p. Thus, the ratio 


e T (I e/p j- 4 2 

e T g/m-p 

is distributed as F(p,m-p), if 0=9 . But from Eqs. 

T ? ° 

(4.18) and (4.19), e g= (m-p)s and it is easily shown that 
e T (I -v(l ) e = (9-e o )X T X(9-0 o ) , if 9 = 0 Q . Hence, 


(0-0 n )X T X(0-0 n ) 

~ ^ F(p,m-p) (4.30a) 

PS 

Examination of Eq. (4.24) shows that the numerator of Eq. (4.26) 
is the difference between the sum of squares regression on 
0 q and the sum of squares of the actual parameters 0. 

High values of F(p,m-p) correspond to rejection of the 
hypothesis that 0 = 0 Q . In particular, to test whether the 

data depends on a specific parameter, 0. of 0, the 

F(p,m-p) ratio is evaluated for 0 Q = 0. 


T 

A sum of squares can be written as y Ay where y is 

vector, A is a matrix of known constants. Then the number 
of degrees of freedom of S. is defined to be the rank of 
A. 1 
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A recursive form of Eq. (4.30a) is 


F(q,m-p) 


(9-0 )X 1 X(0-6 ) 
o' 0 


qs 


(4.30b) 


where q is the number of parameters in the regression at 
any stage in adding (or deleting) parameters to improve fit. 

The F-ratio tests are based on quotients of regression 
"fit" to error "fit." The quantity, 

r 2 = (X9)^(X6) 


= y T y/y T y (4.3i) 

measures the regression sum of squares to observation sum 
of squares. The positive square root of Eq . (4.31), R, is 

the multiple correlation coefficient. The closer R is to 
unity, the better the performance of the subset of regression 
variables. Note that R is the cosine of the angle between 
the data vector and the p dimensional subspace spanned by 
the included subset of regression variables. 

From the sum of squares decomposition (4.20), and 
Eq. (4.31), it may be shown that 

F = ffiP- 1 (4.32) 

(l-R Z )/(m-p) 


2 

which expresses the F-ratio in terms of R . 

The closer R is to unity, the stronger the dependence 
of the data on the regression parameters and the higher the 
F-value. If there is little dependence, R is "small," 
the hypothesis 9- = 0 is "correct," and F is "small." 

This test is dependent on choosing a critical value of F 
which specifies the cutoff levels. This critical F is 
chosen as a function of m (number of observations) and p 
(number of parameters) from tables of F statistical tables 
for a desired confidence level. The confidence level can be 
selected on the basis of a priori knowledge about level of 
noise in data. 
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Significance of individual conditions to the regression 
equations . - The multiple correlation coefficient measures 
the dependence of the data on the complete set of regression 
variables . If several sets of parameters are to be evaluated 
with respect to their respective explanations of the data, 
the R coefficient would serve as one measure of relative 
performance. Unfortunately, it follows from the definition 
of R that increasing the number of explanatory variables 
always increases R, unless the added variables are linearly 
related to other included variables. Beyond a certain point, 
the added variables start fitting the random noise e. For 
example, if the number of variables is increased to m, a 
perfect fit can be obtained and R is one. In order to 
eliminate parameters not significant, only the most linearly 
independent variables are desired. 

Stated in terms of the sum of squares principal, it is 
desired to incorporate those new variables into the regres- 
sion which most reduce the error, given that the other vari- 
ables are included in the regression. If the increase in 
regression sum of squares due to adding a new variable, say 

x. , after variables x- have been included in the regres- 
3 /\ 

sion, is denoted as A| |y| | , then the ratio 


R. 


yx. • x. , x. 


,x j-l ,x j+l» 


A 

X = A l )y 1 1 

’ P e T e + A | | y I I ' 

(4.33) 


measures the partial contribution of 

2 

The square root of R 

yXj a^,..., Ap 


Xj to the regression, 
is the partial correla- 


tion coefficient. Several forms of this definition may be 
written. One is defined by letting y be the residuals 
from the regression of x. on the same variables. It may 
then be shown ^ 


i Ciy *,) 2 

= —1 

(^Kly ) 2 

Note that, in general, ER~~ 1 R 2 . 

j yX j 

Now, let the regression equation y 
titioned as 


(4.34) 


X0 + e be par- 
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(4.35) 


y = X 1 6 1 + X 2 e 2 + e 

where X^ includes q variables and X 2 contains p-q 
variables. Then 


y - X 1 0 1 = X 2 0 2 + e 


(4.36) 


which shows that an estimate of 0^ could be obtained by 

regressing the residuals from the regression of y on X^ 

(i.e., which estimates 0-,). Then the vector y-X,0, is 

i fll 1 1 

regraded as a new observation, say J , which may be 

regressed on X 2 to estimate 0 2 - This decomposition can 
be applied to each possible subset of variables, X^, "bring- 
ing in" new variables from the right to left hand side of 
Eq. (4.36). The requirement on "bringing in" new variables 
may be satisfied by examining the significance of each vari- 
able. 


The F test may be used to determine the significance 

of a single parameter by noting the estimate of the variance 

2 2 2 2 2 2 2 
a , s , is distributed as a X m . p - Hence, s /a ^ (x m _ p ) / (m-p) 

from Eq . (4.19). Then for the parameter 0^, 


0 i -0 i _ (0.-0 i )/ae i 
s 0i s 0i/°0i 


(4.37) 


where Sq^ is the standard error of 0^, which, from Eq. 
(4.20) is 


s 9 i = s y/ ^ii (4.38) 

where s.. is the square root of the ith diagonal term of 

( x T x)-i. 1X 

Since (§ i '0 i )/a i ^ n(0,l), it follows that, by defini- 
tion of Student's t distribution that 
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(4.39) 



t 

m-p 


In particular, if it is desired to test the hypothesis 
0 ^ = 0 (i.e., y does not depend on 0 ), the statistic 

t =■ 0 S 0 ^ is used. It is shown [23] that the F distri- 
bution with 1 and (m-p) degrees of freedom is equivalent 
2 

to the t distribution with m-p degrees of freedom. 

Hence, the significance of individual regression coeffici- 
ents, 0, is determined from F-ratios 

F = 0?/s? . (4.40) 


If the ratio (4.40) indicates a variable is not sig- 
nificant, then the variable is deleted. To bring in another 
variable, the partial correlation coefficients of all other 
parameters are examined. To form the F-ratio for these 
coefficients, Eq. (4.33) (with Eq. (4.38)) may be manipulated 
to show 


R~ 




( 9 : /S 9 i ) 




(m-q) 


(4.41) 


where q is the number of variables already in the regres- 
sion. The corresponding F-test is 


F. 

1 


(m-q) 

-7V- 


(4.42) 


The variable (F-ratio with 1 and m-q degrees of freedom) , 
is calculated for each of the remaining variables. The 
variable with the highest value is then brought into the 
regression. 

Summary . - The variables x^,X 2 »-**»x are postulated 

as possible causative factors in determining the observed 
data, y (Eq. (4.35)). The variables x^,X 2 »...,x are 
ranked according to the magnitude of their individual partial 
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correlation coefficients with y. The variables with the 
highest significance (F to enter) is added to the regres- 
sion. The significance of the added parameter is then 
tested to determine if it is above the critical F value 
for deletion (Eq. (4.40)). As variables are entered, ne\>r 
F-ratios are computed from the remaining variables not in the 
regression since the degrees of freedom have been reduced. 

This process is terminated when no new variables satisfy the 
F- ratio required to enter and when one is to be removed. 

At each variable incorporation, the multiple correlation 
coefficient (R, Eq. (4.31)), the equation F-ratio (Eq. (4.32)), 
and the standard error of the entered parameter (Eq.(4.38)) 
are all included. 

It is noted that, in general, the particular variables 
finally selected are not unique. Use of orthogonal variables 
would result in uniqueness. 

Numerical methods . - The previous sections have demon - 
strated that both classical and subset regression parameters 
are obtainable from steps in the solution of a set of linear 
equations (ref. Eq. (4.15)). In order to reinforce this 
connection, consider the augmented data-coeff icient matrix 
(with rank p+1) 


A 


T T 

X*X X x y 


y T X y T y 


Using the inversion of Eq. (4.15), it is found that the 
- 1 f 

final element of A , a , is 


a £ - 


y T y - y T X(X T X)' 1 X T y 


T T T 

y y - e^X 


(1 - R 2 ) y T y 

i.e., the multiple correlation coefficient (and the sum of 
squares | | e | | ) are directly obtainable. Similarly, it is 

found that the diagonal element of (X T X)’ 1 corresponding 
to the ith coefficient of x, is 
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1 




Cl - P 2 )x T x. 

1 1 


2 

where is the multiple correlation coefficient of the 

regression of on Xj ,x 2 , . . . ,x i _ 1 ,x i+1 , . . . ,x . It may- 

be further shown that the variance of 0^ is 


/v y t — "i 

var 0. = a^(X i X) i ^ = 


(l-P?)x T x. 


Computationally, the inversion of A is based on the 
Gauss-Jordan pivot algorithm. Let the kth diagonal element 
of A be nonzero. Applying a Gauss- JorcTan pivot on this 
element is a new amtrix whose ijth element is a.., i.e., 


! a ij " a ik a kj /a kk 
' a ik /a kk 
a kj /a kk 
1/a kk 


i f k, j f k 
i + k, j = k 
i = k, j f k 
i = k, j = k 


The final result of this inversion is the matrix B, 


B 


(X T X) " 1 0 



The recursive algorithm of the Gauss-Jordan pivot sweeps 
through the A matrix, generating statistical parameters 
which guide the deletion and addition of new variables. 

Details of this selection for the Gauss-Jordan pivot are 
found in Ref. 24. One further parameter of this method is 
the tolerance. This is a parameter of the computation algo- 
rithm itself and is a measure of the accuracy of the calcu- 
lations and computer. If a^ is the value of the ith 

diagonal element of the non-negative matrix ) after 

pivoting on several elements, then a ij/ a ii as t ^ ie tolerance. 
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CHAPTER V 


ROTORCRAFT PARAMETER ESTIMATION 


Introduction 

Rotorcraft parameter identification is the process of 
extracting the stability and control derivatives for a given 
mathematical model structure using a set of flight test or 
simulated data. In the past, aircraft identification has 
been carried out as an ancillary part of experimental flight 
test work. The principal source of stability and control 
derivatives has always been wind tunnel results and theoreti- 
cal calculations. However, as the result of what is often 
gross disagreement between wind tunnel and flight test 
derivatives, and the known difficulties of obtaining wind 
tunnel values, as well as extrapolating these values to full 
scale, there is an increasing emphasis on obtaining these 
parameters from actual flight data. 

In the field of system identification, most of the ef- 
fort had been, for a number of years, confined to the solu- 
tion of low order linear problems. However, during the past 
few years, powerful digital techniques have been developed 
which are capable of solving the more difficult nonlinear 
aircraft problems. The problem of extracting aerodynamic 
coefficients for rotorcraft is one of the most difficult 
identification problems. 

The development of accurate parameter estimation is a 
logical step following model structure determination. This 
is because the model structure estimation step (described in 
the previous chapter) gives a reasonable model of rotorcraft 
forces, moments and other dynamic parameters, but the esti- 
mates of parameters resulting from the least square approach 
may be of questionable accuracy. Several procedures have 
been used for estimating nonlinear system parameter in the 
past. Because of the complexity of rotorcraft dynamics, the 
selection of an appropriate identification method is extremely 
important. This chapter discusses various identification 
methods pointing out the advantages and disadvantages of 
each. Finally, one method is selected for further investi- 
gation. 


Rotorcraft Parameter Estimation Approaches 

General discussion .- Most identification methods that 
have been or are currently being applied to rotorcraft prob- 
lems can be classified as either: 
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(1) equation error methods, 

(2) output error methods, or 

(3) advanced methods. 

These methods differ by: (1) the performance criterion that 

they are developed from, (2) the kind of estimates they pro- 
duce, and (3) the problems to which they can be applied. 

Equation error methods [25] assume a performance cri- 
terion that minimizes the equation error squared (the process 
noise). All of these methods are basically least squares 
techniques. In general, it is necessary to measure the 
response variables and their time derivatives to apply these 
methods. This results in n or more linear equations in 
n unknowns. For the case where the time derivatives are 
not measured, various "method functions" are used to operate 
on the equations (e.g., time derivatives, Laplace or Fourier 
transforms, etc.), to obtain equations that are linear in the 
unknowns. Since these methods do not allow for measurement 
errors (sensor noise, etc.) they result in biased estimates 
when this type of error does exist. The principal use of 
these methods are as start-up techniques for the output error 
and advanced methods. The equation error methods have been 
used or are being used by Calspan Corp. [25], Air Force 
Flight Test Center [26], and by Gerlach [27]. 

Output error methods [28-33] minimize the square of the 
error between" the actual system output and the output of a 
model. These methods assume measurement noise but no process 
noise. Typical output error methods include Newton-Raphson; 
Gradient methods; the Kalman filter (without process noise); 
and modified Newton-Raphson, differential correction, and 
quasilinearization (all three of which are the same method). 

It should be noted that analog matching is basically 
a manual form of the output error methods (it attempts to 
match the simulated response to the actual response) . This 
method has been used by the Air Force Flight Test Center [26], 
the Naval Air Test Center [34], and the NASA-Dryden Flight 
Research Center [35] for the F-104, X-15, B-70, HL-10, M2/Fe, 
X-24, and PA-30 aircraft. The modified Newton-Raphson method 
has been used extensively in flight test applications for the 
past several years. It is the one method that has been used 
on an operational basis and for which the most experience exists. 
It is known that this method has been or is being used by: 

(1) the NASA Dryden Flight Research Center [35] on lifting 
bodies (X-24, F-3, X-14, SB-70, 990, HL-10, M2/F3) , Jet Star 
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and PA-30 aircraft, (2) the NASA Ames Research Center [36] 
on the Learjet, XV- 5, 990, and C-8 aircraft, and (3) the 
NASA Langley Research Center [37] on the XC-142, Navion and 
F-4 aircraft. 

The principal disadvantage of output error methods is 
that, because they do not include process noise in their 
performance criterion, the results degrade when process 
noise (gusts, modeling errors) exists. This may result in 
the program not converging or in estimates that have large 
variances or poor estimates [38], However, as long as these 
methods are applied to linear flight regions, or where the 
form of the equations are known, and where gusts are not 
significant, they work very well. 

There are two advanced methods which have been applied 
to flight data with both process and measurement noise. 

These are: 

(1) The Kalman filter and smoother [39-41] 

(2) The maximum likelihood method [42-47] 

Because of the advanced nature of these algorithms, they are 
discussed in some detail in the following subsections. 

Extended Kalman filter/smoother parameter identification 
methods . - In 1961-1962, R. E. Kalman published a sequence of 
papers in which he proposed a filter with the capability of 
estimating randomly disturbed states of a system from noisy 
measurements. The Kalman filter has, since then, been the 
object of intense research and development. It is signifi- 
cant to note that the Kalman filter for state estimation is 
a particular solution to a more general problem of estimating 
both unknown states and parameters of a dynamical system. 
Kalman, himself, suggested in this early work that the fil- 
ter could be adapted for estimation of the parameters as 
well as the states by simply augmenting the parameters to 
the state and estimating both simultaneously. This is still 
the first approach that is usually tried in developing para- 
meter identification algorithms, because it does have a 
strong intuitive and simple implementation. Unfortunately, 
the approach is found to suffer from certain inherent prob- 
lems in actual implementation, particularly for high order 
systems perturbed by both process and measurement noise. 

There is a simple reason for the difficulties encounter- 
ed by this approach. The Kalman filter is strictly a state 
estimator derived on the assumption of known parameters. 
Assuming a filter of the Kalman type for also estimating the 
parameters is an "after-the-fact" adjustment. Practically, 
the result is that to apply the Kalman filter to even a 
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linear system requires that the linear system be converted 
into a nonlinear system. Although this transformation does 
not, in itself, degrade the algorithm effectiveness, it does 
cause other problems which can. Other problems which are 
reported are; 

(1) The Kalman filter approach, even with the algo- 
rithm modifications such as smoothing, produces 
biased estimates of the parameters [41]. Hence, 
the parameter estimates may be in error. 

(2) Any Kalman filter algorithm requires initial 
guesses for both states and parameters, as well as 
their respective state and parameter variances 
[39,41]. The necessity for providing good state 
variances is not a severe limitation. However, 
parameter variances are usually much harder to 
guess than state variances. Unfortunately, the 
Kalman filter approach does require accurate para- 
meter variance initiation to obtain good estimates. 
A common technique is to use initial guesses and 
variance obtained from a least squares regression 
on the data. These results are biased, however, 
and in order to account for bad guesses, arbitrary 
variance multiplication "factors" must be used to 
ensure even algorithm convergence [41] . Such 
"guesses" for each type of case run require highly 
skilled personnel. 

(3) The Kalman filter approach is often claimed to 
estimate parameters in the presence of process 
noise, even though there is no theoretical justi- 
fication for its ability to do so. What is done 
is actually to either filter the data by some 
smoothing technique or simply adjust the process 
noise covariance terms of the filter until con- 
vergence is achieved. Molusis [39] uses the 
smoothing approach to filter the data. The prob- 
lem with this approach is that the resulting data 
may contain less information about the parameters 
of interest, again causing errors in the estimates. 
The iterating smoothing approach suggested by 
Chen, et al . [41] is probably the best implementa- 
tion of extended Kalman filter for parameter 
estimation. 



Maximum likelihood methods .- The previous section 
alluded to the fact that the Kalman filter for state estima- 
tion is a particular solution of the more general problem of 
simultaneous state and parameter estimation. Realizing the 
difficulties of trying to use the Kalman filter approach in 
practical application, a return to the original formulation 
of the simultaneous state and parameter estimation problem 
leads to a more general estimator. This formulation- -a 
maximum likelihood approach- -is one of the ways a Kalman 
filter state estimator is derived. The maximum likelihood 
method is illustrated in Figure 5.1. Conceptually, this tech- 
nique can be summarized as follows : 

"Find the probability density functions of 
the observations for all possible combina- 
tions of unknown parameter values. Select 
the density function whose value is highest 
among all density functions at the measured 
values of the observations. The correspond- 
ing parameter values are the maximum likeli- 
hood estimates." 

The method is a combination of three steps: 

Cl) A Kalman filter to estimate the states. 

(2) A Levenberg-Marquardt optimization technique 
for the parameter estimates. 

(3) An algorithm to estimate the noise statistics 
(mean and variance of the measurement and process 
noise) . 

The method has the following important features: 

(1) Parameter Estimate Variances . In estimating para- 
meters from flight data, it is essential to pre- 
dict the error in these estimates, that is, the 
confidence level associated with the estimates. 

The maximum likelihood method automatically com- 
putes the lower bounds on the variances of the 
parameter estimates. Since the maximum likelihood 
estimate is asymptotically efficient, the actual 
variances in the parameter estimates approach 
these bounds for long data records. The method 
does not require initial parameter variances to 
initialize the algorithm. 
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FLIGHT TEST DATA, WIND TUNNEL VALUES OF 
AERODYNAMIC PARAMETERS 



FIGURE 5.1.- FLOWCHART OF MAXIMUM LIKELIHOOD 
IDENTIFICATION PROGRAM 










(2) Noise Statistics . Many identification algorithms 
assume that the means and variances of the process 
and measurement noise are known. In fact, these 
quantities are generally not known and should be 
identified. Typical noise parameters that should 
be identified include sensor biases, c.g. off- 
sets, the variance of the random measurement noise, 
and the mean and variance of the 

require a priori knowledge of the process or meas- 
urement noise covariances but determines them as 
part of the identification procedure. 

(3) State Estimates . In addition to estimating the 
parameters of the aircraft models (stability and 
control derivatives) , it is also necessary to 
estimate the airplane response variables (state 
variables). The maximum likelihood method obtains 
the best mean square state estimates as part of 
the parameter identification process. 

In addition, the generality of the maximum likelihood method 
has important operational features: 

(1) In the absence of process noise, the method is 
equivalent to the output error method which many 
flight test agencies already use. 

(2) In the absence of measurement noise, the maximum 
likelihood method reduces to a least squares 
approach. 

(3) The method requires no "fine tuning" or "fiddling" 
to provide estimates for new cases. 


Rotorcraft Parameter Identification Method 


It is clear from the discussions of the previous sec- 
tion that the maximum likelihood provides an optimal blend 
of optimality and flexibility in rotorcraft identification 
problems. It is for this reason that the maximum likelihood 
method has been selected for investigation in rotorcraft 
applications. The maximum likelihood method is one of the 
most flexible techniques in statistics for identification of 
parameters from input-output data. Suppose it is possible 
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to make a set of observations on a system, whose model has 
p unknown parameters 6. For any given set of values of 
the parameters 0 from the feasible set 0, we can assign 
a probability p(Z|0) to each outcome Z. If the outcome 
of an actual experiment is z, it is of interest to know 
which sets of values of 0 might have led to these observa- 
tions. This concept is embedded in the likelihood function 
SBkQ\z). This function is of fundamental importance m 
estimation theory because of the likelihood principle of 
Fisher and others [48-50] which states that if the system 
model is correct, all information about unknown parameters 
is contained in the likelihood function. ^The maximum like- 
lihood method finds a set of parameters 0 to maximize this 
likelihood function 


0 = max SB [0 | z) (5.1) 

6e0 


In other words, the probability of the outcome of z is 
higher with parameters § in the model than with any other 
values of parameters from the feasible set. Usually it is 
more convenient to work with the logarithm of the likelihood 
function (it is possible to do so because the logarithm is 
a strictly monotonic function) . 

Consider a nonlinear rotorcraft of the following form 


x = f(x,u,9,t) + r(0,t)w 0 £ t < T 

E[x(0)] = x 0 (e) 

E{[x(0) -x o (0)][x(O) -x o (0)] T } = P o (e) 


where 

x(t) is nxl state vector 
u(t) is £xl input vector 
0 is pxl vector of unknown parameters 

T is nxq process noise distribution matrix 

w(t) is qxl random process noise vector 
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taken at discrete times 


Let there be m measurements y(t, ) 

t K 

r k 

y(t k ) = h(x(t k ) ,u(t k ) ,0,t k ) + v(t k ) 

k = 1,2,3,. . . ,N 

w(t) and v(t k ) are Gaussian random noises with the follow- 
ing properties 

E[w(t)] = 0 E[v(t k )] = 0 E[w(t) v T (t k )] = 0 
E[w(t)w T (x)] = Q (9 , t) <5 (t-r) , 

E[v(t :j )v T (t k )] = R(9 ,t^. ) « j k 

The log- likelihood function for this mathematical model can 
be shown to be of the form 


log[<2?(e|z)] = 4 Z (v T (i) B~ 1 (i) v(i) + log I B CiD ! } 
Z i = l 


where v(i) is the innovations vector at sample point t^ 

and B(i) is its covariance. An estimate of the unknown 
parameters is obtained by maximizing the likelihood function 
or the log- likelihood function from the feasible set of para- 
meter values 


0 = max log [ c ^(0 | z) ] 

0E0 


(5.2) 


= max 


N 

\ Z {v T (i)B _1 (i)v(i) 
L i = l 


+ log | B (i) j } 


(5.3) 


The log- likelihood function depends on the innovations 
and their covariance. To optimize the likelihood function, a 
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way must be found for determining these quantities. Both 
innovations and their covariance are outputs of an extended 
Kalman filter. The developments of the required Kalman 
filter is discussed in subsequent paragraphs. 

The extended Kalman filter is conventionally divided 
into two parts. In the first part, called the prediction 
equations, the state equations and state estimate covariance 
equations are propagated in time from one measurement point 
to the next. In the second part, called the measurement 
update equations, the measurements and associated measurement 
noise covariances are used to improve state estimates. The 
covariance matrix is also updated at this point to reflect 
the additional information obtained from the measurements. 

Prediction equations .— The state prediction is done using 
the equations of motion. Starting at time t^_^ with current 

estimate x(i-l|i-l) of the state x(t^) and the covariance 
P(i-l|i-l), the following equations are used to find the pre- 
dicted state x(i|i-l) and the associated covariance P(ifi-l); 
see Bryson and Ho (Ref. 10). 

^ x(t|t i _ 1 ) = f (x(t |t i _ 1 ) ,u(t) ,0,t) (5.4) 

Pftjt.^) = F (t) p(t|t i _ 1 ) + P(t|t i _ 1 ) F T (t) + rQr T (5.5) 


The nxn matrix F is obtained by linearizing f about the 
best current estimate 


F (t ) 


3f (x(t |t i _ 1 ) ,u(t) ,0 ,t) 
3x(t |t i _ 1 ) 


(5.6) 


Using (5.4) to (5.6), we can obtain 


x(ti|ti_i) A x (i 

1 i- 1) 

(5.7) 

p (t i |t i _i) A P(i| 

|i-l) 

(5.8) 


Thereafter, the measurement update equations are used. 



Measurement update equations .— The covariance and state 
estimate are updated using the measurements. The necessary 
relations are derived by Bryson and Ho (ref. 10) and are 
presented here without proof. The innovation and its covari- 
ance are 


v(i) = y (i) - h (x(i | i-1) , u(t i ), 9, t ± ) (5.9) 

B (i) = H(i) P ( i | i - 1) H T (i) + R (5.10) 

where H is obtained by linearizing h, 

8h (x(i | i-1) ,u(t . ,6 ,t. )) 

H(i) = — (5.11) 

3x(i | i-1) 

The Kalman gain and the state update equations are 

K(i) = P(i | i-1) H T ( i) B" 1 (i) (5.12) 

x(i.|i) = x(i|i-l) + K(i) v(i) (5.13) 

Finally, P(i|i), the covariance of error in updated state, is 
obtained by 


P(i I i) = (i-K(i) H(i)) P (i | i-1) (5.14) 

The maximum likelihood method can be shown to reduce to 
output error in the absence of process noise. If there is no 
measurement noise, it reduces to the equation error method. 

These simplifications are discussed in Appendix A. 

Optimization procedure .— It is possible to use any one 
of several numerical procedures for this optimization problem. 
The authors have found by experience that the Levenberg- 
Marquardt optimization technique provides the most rapid con- 
vergence. This technique has both the Newton-Raphson and the 
steepest descent algorithms as special cases. The general 
Levenberg-Marquardt approach is discussed in references [51 
and 52], The algorithm requires computation of the first and 
second order partials of the log- 1 ikelihood function. 
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The last three terms in the equation for second partial of the 
the log- likelihood function involve second partials of innova- 
tion and its covariance. Those terms are usually dropped. 

So, the second partial is approximated by 


3 Z log(g(9 lz)) _ 


N 


30 j 30], 


(i) B " 1 m 

1 J 30 k 


T , . 

\> (l 
30v 


z 1 

i=r 
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CD ' 
CD 
t-J. 
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(i) 

3B(i) 

30, 


- B - 1 (i) B _1 (i) v (i) 


30 


30i 


(5.16) 


+ v T (i)B _1 (i) B _1 (i)v(i) 
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B-l(i) Mill B~ 1 (i) MAi 
^ J 30 j 30 k 


The gradients of innovation and its covariance for parameter 

0 . are : 

3 


3v(i) _ 


30 


3 


3h 

3x 


x=x(i | i-1) 


3x(i | i-1) _ 3h_ 
30 " 30 


(5.17) 


j 1>2,... ,p; i 1»2,...,N 


p Ci I i-l)H T Ci) + H(i) - - P( y i ~- 1) H T (i) 

d0j 3 0 j 30 j 

. H( i) P(ili-i) 2!&±1 + JB- 


(5.18) 


30. 


30. 


3 i 
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Recursive equations can be obtained for gradients of the 
predicted state and its covariance. This is done in stages 
by using the prediction and measurement update equations of 
Kalman filter. Differentiating (5.4) to (5.6) with respect 
to 0 . 


d 9x(t|t i _ 1 ) 9f (x(t |t i _ 1 ) ,u(t) ,t,e) 


art 


99 


90. 


3 WV J 

9f (x(t|t i _ 1 ) ,u(t) , t , 0 ) 9x(t|t i _ 1 ) 


9x(t | t i _ 1 ) 

9x (0 1 0) _ ax o( 9) 
90 j 90 j 


99 j 


A. 9P(t|t i-l } . 3F(t) p(t , t ) + F(t) 

dt 30- 39- * r i-l J 


3P(t|t i _ 1 ) 

90^ 


(5.19) 




+ r -AL. r T + tq 

90 j W 90J 

9P(0 | 0) _ 9P o^ 


(5.20) 


90. 


90. 


Vl ^ 1 * 


(5.21) 


9F(t) 

90 


9 2 f (x(t |t i _ 1 ) ,u(t) ,0,t) 

9 9 j 3xTTlt~7) 


j l»2,...,p 


The sensitivity functions are updated at measurement points 
by differentiating (5.11) to (5.14) with respect to 0^ . 
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(5.22) 


3H(i) _ 3 Ci | i - 1) >u(t i ) ,0,t i ) 


30- 


36^ 3x(i|i-l) 


IMil . H T (i) B-Vi) * P(ili-l) 4^ 

d y j d y j d u j 


- P(i|i-1) H T (i) B -1 (i) B~ 1 Ci) 

aej 


(5.23) 


3x(i I i) = 3x ( i | i- 1) + 3K(j) + 3v(i) 


30- 


30 j 


30 j 


30 


(5.24) 


9P ^ |i} ■ I-K(i) H(i) - P( ^ 1 ~ 1 - ) - - H ( i) P(i|i-1) 


K(i) 


(5.25) 


j 1*2,3, ...,p 


The negative of the matrix of second partials of the log- 
likelihood function is called the information matrix M. The 
step size A0 for parameter estimates is given by 


A0 


p(M+aI) -1 8 lo g ( g eqM 

3 0 


(5.26) 


where a is the Marquardt parameter and p is a scaling term. 
Observe that as a -*■ o the step approaches that of the Gauss- 
Newton procedure. As a -* ® (with p/a finite) the step 
approaches that for the steepest descent algorithm. During 
the optimization, it is desirable to change the Marquardt 
parameter to promote rapid convergence. If it is found that 
the step is too large (i.e., the log- likelihood cost increases), 
the step is cut by increasing the Marquardt parameter by a 
given factor and another step is taken. If the cost decreases, 
the Marquardt parameter is decreased by the given factor so 
that the following step is larger. This is to prevent the 
slow convergence that results from taking a small step — as 
is often characteristic of steepest gradient algorithms. 
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Linear systems .— In a linear system, the functions f 
and h are defined as 


f(x,u,0,t) F(8,t)x + G(0,t)u 

and h(x,u,0,t) _A H(8,t)x + D(0,t)u 

The basic algorithm is the same. However, some of the 
equations can now be simplified. Equations (5.4) and (5.9) 
now become 

^■x(t|t i _ 1 ) = F (0 , t) x(t |t._ 1 ) + G(0,t) u (t ) (5.27) 

and 

v(i) = y (i) - H(0,t i ) x (i | i- 1) - D(0,t.) u(tj) (5.28) 


Equations (5.17) and (5.19) can be written as 

lii. . H (8,t.) aiciiiii . idii-D 

0-: 904 80-; 1 


3v 

90 


J 


9D(0,t.) 


C5.29) 


and 


axftlt.^) _ 9x(t|t i _ 1 ) 9FCe>t) c 

F(e,t) —90— x tt|t._ 1 ) 


dt 


+ 


9G (8 ,t) 
90 j 


u(t) 


(5.30) 


All other equations remain the same. 


There is considerable reduction in computation require- 
ment for time- invar iant linear system. In this case, matrices 
F, G, H, D, T, Q and R and their derivatives with respect to 
parameters are constant. 
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For many cases, the Kalman filter is in steady state for 
the duration of the experiment. This occurs when the Kalman 
filter is in operation for a sufficiently long time and the 
process and measurement noise covariances do not change. The 
Kalman gain and the innovations and the state covariances ap- 
proach constant values. The time update and measurement up- 
date equations for the covariances are 

( 

57 P(t|t i _ 1 ) = FP(t|t i _ 1 ) + P(t|t._ 1 )F T + rQr T (5.31) 


K = P ( i i i- 1) H 

T - 1 
B 1 

(5.32) 

B = HP (i | i- 1) H T + R 

(5.33) 

P(i | i) - (I-KH) 

P(i | i-1) 

(5.34) 

By definition of the steady state 

P(i-1 1 i-1) = 

P(i 1 i) 


Therefore, from (5.31) 

P(i | i-1) = e FAt P(i- 1 | i-1) e F At + 1 

f h TQT T f^h, 


J 

t 

_T 

:i-l 

ft-i . _ _T. 



= e FAt (I-KH) P(i|i-l)e F At + J e F( ' ti “ T ' ) rQr T e F ^"^dx 


A <KAt) (I-KH) P ( i J i-1) 4> (At) + Q" 

= <j> (At) (P(iji-l) - KBK T ) (j) (At) + Q" 


(5.35) 


Using (5.32), (5.33) and (5.34), we can solve for P(i|i-1) 
and then find K and B. Also, it can be shown that 


3 P 
30 


= A-, 


3 P 
36 


A^ + A 2 " <h PA^ P0 


(5.36) 
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3K 


(5.32) 


= (I-KH) h T (HPH T + r)" 1 + A 2 

= .liL ph T + H H T + HP 
30 • 30 j 30j 30 j 30 ■ 

where 

A x = <j> (I-KH) 

A 7 = -ii- (I-KH) P<j) T + <|) (I-KH) P <j)K — Pcf> T + ^21 

2 30^ J ’ 30 j y 30 j v 30^ 

A = HL b _ 1 H + H T B _1 — PH T + HP |S- + JiL B _1 H 

3 38 j 3 0 j 30 j 30 j 

P A P(i | i-1) 


(5.38) 


Thus, it is possible to solve for 


3P 

90. 


using (5.36) and then 


find and from (5.37) and (5.38). Equation (5.36) is 

30 j 


30 . 

3 

linear equation in 


3P 

3 e. 


and the coefficient of the unknown 


matrix does not depend on the parameter 0^ . Thus, the 


sensitivity of state covariance matrix can be determined 
very quickly for all parameters. Once the sensitivity of P, 
K and B for unknown parameters is determined, only state 
sensitivity equations need to be updated. 


a 


An approximation simplifies the problem further. The 
unknown parameters are defined to include elements in K and 
B matrices instead of Q and R. Optimizing the log-likelihood 
function for parameters in B gives 


B 



N 

Z v (i) 
i=l 


v (i) 


(5.39) 


1 

1 

J 

1 



1 


n 


1 


i 
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t 
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The gradient of the log-likelihood function with respect to 
other unknown parameters is 


3 Log (g(6 I z) ) _ 
36 j 


I v T (i) B' 1 Mil 
i-1 30 j 


(5.40) 


The sensitivity of innovations to parameters is determined 
using the following recursive equations 


d ~ 


— x(t|t i _ 1 ) = Fx(t|t i _ 1 ) + Gu (t) 


(5.41) 


d 3F - 


dt 


36 


36 


J J 

3 = 1 > 2 , . . . , p 


- x(t|t i _ 1 ) + F 




39 j 


36j 


u(t) 


t. < t < t. 
l-l l 


(5.42) 


9v(i) 

30 4 


v(i) = y (i) - Hx ( i | i - 1 ) - Du(t.) 


x (i | i) = x ( i | i-1) + Kv (i) 


3 0 j 1 39 j 


SxCij i) = 3x(il i-1) _3K_ m K 3v(i) 

304 304 304 304 


) 

(5.43) 


(5.44) 

Be, “"i 5 

(5.45) 

3v(i) 

304 

(5.46) 


Note that 
3K 


30 


3 


3 1,2, ...,p 


= 0 if 9 . is not an element of K matrix 
3 


- if 9 j A K.^ 


where 1^ is a matrix of all zeroes except a 1 at the j' ,k' 
pos it ion . 
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CHAPTER VI 
INPUT DESIGN 


In previous sections, we have considered aspects of rotorcraft 
system identification methodology which address the post flight 
data processing phase. Data are presented to the procedure, and 
some level of model and coefficient estimation achieved. One does 
not have to perform this task many times before it is evident that 
no post-flight data procedure can identify coefficients if the in- 
formation on these parameters is not in even the uncorrupted mea- 
surement response. An input design requirement, therefore, has 
arisen because of the need to improve the efficiency of flight 
testing by obtaining more accurate estimates from response data in 
less time. 


Requirements for Input Design 

There are several factors that must be considered when choosing 
inputs for flight tests. These include: 

(1) Efficiency: The inputs should improve the efficiency of 

stability and control flight testing by maximizing coefficient 
identif iability in fewer maneuvers. 

(2) Structural and Safety Constraints: The inputs should ex- 

cite the modes of interest without violating rotorcraft safety and 
structural constraints. 

(3) Identif iability : The inputs selected must improve param- 
eter estimation accuracy and avoid identifiability problems in the 
post-flight data analysis stages. 

(4) Pilot Acceptability: If the flight test is to be carried 

out with a pilot onboard the rotorcraft, it is necessary that the 
control inputs be acceptable to the pilot. The inputs should not 
maneuver the aircraft into a flight regime from which a pilot can- 
not recover. In addition, the inputs should be reproducible by the 
pilot . 

(5) Instrumentation: The inputs must consider specific in- 

struments available, and their dynamic range and accuracy. The 
primary impact of the instruments on input design is on the signal/ 
noise ratios which the response must have for sufficiently accurate 
data . 


(6) Modeling Assumptions: The inputs that are designed must 

also consider the model that is assumed. For example, inputs 
chosen for a linear mode should not cause such large rotorcraft 
motions that the assumption of constant stability and control de- 
rivatives is invalid. 
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Input Design Methods 


Since the first efforts of applying parameter extraction 
technology to aircraft flight test data, many different control in- 
puts have been used. In one approach, the flight vehicle is ex- 
cited by sinusoidal inputs over a range of frequencies, usually 
around the natural frequency of the specific mode, until a steady 
state is reached at each frequency. The parameters of a suitable 
linear model are selected to obtain the best fit to the variation 
with input frequency of the output/input amplitude ratio and phase 
difference. These inputs work satisfactorily but require much 
flight test time. Pulse inputs have been used to identify simple 
low order linear systems. Doublets, steps and finite duration 
pulse inputs have also been used to identify aircraft parameters. 
However, the estimates of certain parameters may be quite poor and, 
in some cases, a set of parameters may not be identifiable at all. 

The definition of information matrix, M, which provides a 
quantitative meaning to the knowledge about a certain set of param- 
eters, is the starting point for a systematic study of the input 
design problem. The inverse of the information matrix, referred to 
as dispersion matrix, is a bound on parameter error covariances, 
i . e . , 

cov j (6 -0) (9 -0) T J >_ M _1 A D 

A 

where 0 is the actual value of a parameter and 0 is the esti- 
mate. Most of the analytical methods in input design use a func- 
tion of information or dispersion matrix as the extremizing crite- 
rion . 


Optimal input design for dynamic systems .- In this section, we 
describe a method consisting of two different techniques for design 
of input signals which provides estimates of unknown parameters in 
linear systems. The first technique used the time domain represen- 
tation of system dynamics and develops methods to compute the time 
history of control input sequence for any duration of the experi- 
ment. In the second technique, the system is assumed to be in 
oscillatory steady state and a frequency domain representation is 
utilized. This gives the optimal control input spectrum. The 
corresponding time history is "optimal" only for long experiments. 
The time domain approach is computationally much more complicated 
than the frequency domain approach. The two approaches are, there- 
fore, complementary. The frequency domain approach is suitable for 
long experiments and the time domain approach should be used for 
short and medium duration experiments. 

Consider a linear, time - invariant , dynamic system described by 
the following differential equation: 
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X 


Ax + Bu 


( 6 . 1 ) 


x(0) =0 0 £ t <_ 1 

where x is a nxl state vector, u is a qxl input vector, and 
A and B are appropriate matrices which depend on m unknown 
parameters 0. Let there be noisy measurements of p linear combi- 
nations of the state vector 


y = Cx + v 


( 6 . 2 ) 


where v is assumed to be a Gaussian white noise vector with in- 
tensity R, and matrix C is a function of 0. Since we are 
mostly concerned here with finite time problems, we may assume 
without loss of generality that the final time is one. For prac- 
tical reasons, we impose a constraint on a quadratic function of 
the input and the states 


J = ( (x T Qx + U T U) dt < E (6.3) 

Jo 

For linear sytems , the inequality sign can be replaced by a equal- 
ity sign. The information matrix, M, for parameters 0 is 



The input u(t) is chosen to minimize parameter estimation errors 
(in terms of the weighted trace or the determinant of the disper- 
sion matrix) subject to the constraint (6.3). This problem can be 
reformulated [53] as a two-point boundary value problem. 


d 

X 


A -y B B T 


X 

dt 

X 


C T R _1 C -A T 


X 



(6.5) 


x (0) =0 , A (1) =0 
~T 

u opt = -yB X 


( 6 . 6 ) 


The smallest value of constant y for a nontrivial solution to 
the two-point boundary value problem gives the optimal control 

input. The quantities x, A, B, T and R are defined below: 
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R 0 . . . 0 

R = o R 

• • • 

• • • 

o • • • R 


( 6 . 11 ) 



The solution to this two-point boundary value problem, based on the 
symplectic properties of the Hamiltonian matrix H, is discussed in 
[ 54 ]. 

If the duration of the experiment is much longer than the sys- 
tem characteristic time, it is possible to design inputs based on a 
variety of criteria quickly by making the assumption of steady 
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state. In an ingenious approach, Mehra [55] converts a linear time- 
invariant system into its frequency domain representation. This 
eliminates the dynamics of the system. The parameter estimation 
problem becomes a regression problem in which input frequencies 
and the power in each frequency are the control parameters. These 
parameters, which define the input design, are chosen by an itera- 
tive procedure. 


Consider the state space representation of a discrete -time 
linear system 


x(k+l) = $x(k) + Gu(k) k=l , 2 , 3 , . . . ,N 
and the noisy measurements 

y (k) = Cx (k) + v(k) 


( 6 . 12 ) 


(6.13) 


$, C, and H are appropriate matrices and contain m unknown 
parameters 0. Fourier transform y(k) to get 


-i Ijrn , „ 

y(n) = C(e J N - $) b u(n) + v(n) 
AW (n , 9 ) u (n) + v(n) 


(6.14) 


As the number of sample points increases, the information matrix 
per sample approaches [55] 


„ _ 1 Dq f 9W* c -l r . 3W , v r , 

M ~ Tn Re \ 30 S w M 5? dF uu (w) 

J-l T 


(6.15) 


where F is the spectral distribution function of u and S 

uu r vv 

is the spectral density of v and superscript denotes the con- 

jugate transpose of a matrix. An algorithm based on this approach 
is outlined in [54]. 


Optimal multistep inputs for parameter identification .— Several 
algorithms have been developed for a solution of the boundary value 
problem leading to the solution of the optimal input design problem. 
The application of this theory has been limited to low order time 
invariant systems due to the complexity of the computation proce- 
dure. The optimal input selection algorithms may be considerably 
simplified, for time variant as well as time invariant systems, if 
the inputs are restricted to be multistep. In this approach, 
the state equations are first propagated to obtain the state time 
histories. A simplifying assumption is then made by approximat- 
ing the states to be multistep and are represented as sums of 
Walsh basis functions. They may then be determined by using the 
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methods developed by Corrington [56] and applied to the linear quad- 
ratic feedback control design problem by Chen and Hsiao [57] . This 
approximation may degrade the quality of the input, particularly 
if the state variables change significantly over the duration of a 
step. 


We consider the system defined by equations (6.1) and (6.2) 
subject to the constraint (6.3). Matrices A, B and C may be 
functions of time. 

The information matrix for parameters 0 is 

0[C(t,9) x ( t) ] ) T R -1 OrC(t,9) x (t)] 

| 99 ) ( 90 

the augmented vector x are defined by the equation 

x + B U (6.17) 

are defined by equations (6 . 8) - (6 . 10) . Further, 

= x(t) + C(t,0) . (6.18) 

i i 

A typical element of the information matrix may be written in terms 
of x(t) 

1 T 

M Aj = x(t) TT(t)R' 1 T j (t)x(t)'dt (6.19) 

where T^ and T^ are defined by equation (6.10). 

It is shown that certain mild restrictions on the class of in- 
puts simplifies the input design problem considerably. Two cases 
are considered in this report. 

Case 1 : The inputs are restricted to be multistep with s 

steps each of duration A. The state and output sensitiv- 
ities are obtained by a solution of the differential equa- 
tions . 

Case 2 : The inputs are restricted to be multistep with s 

steps each of duration A. The state and outputs are 
approximated to be multistep. 

As the number of steps is increased the multistep inputs approach 
the general unrestricted input. 


j dt (6.16) 


M =f 1 

Jo 

The dynamics of 
X = A 

where A, B, x 

3z (t) 
39 , 
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The equations for the information matrix and the quadratic 
constraint are specialized for multistep inputs. These will be 
used in the derivation of the multi-step input algorithms. The 
two cases are discussed separately because different techniques 
are involved in determining the augmented state vector x(t) . 


Case 1 


Under the first approximation, a general representation for 
input u is 

u(t) = Hs (t) (6.20) 

where H is a qxs matrix and s is an sxl vector such that 
s^(t) = 1, (i-l)A<_t < i 
= 0, elsewhere. 


( 6 . 21 ) 


Let Xj^ (t) be the response of the augmental sensitivity system 
(6.17) when a single element H — of H equals one and the re- 
maining elements are zero. In the general time-varying case x-j^ (t) 

may be obtained by a direct solution of the differential equation 
for each i and j . Then because of the linearity of the sensi- 
tivity system, it is clear that 

q s 

x(t) = E Z H. . x. . (t) 
i=l j=l 1J ^ 


( 6 . 22 ) 


qs 

E h. x 1 ? (t) 

i=l l i y J 


where 


(j -i)q+i ^ ^ i j 


Using Equations (6.19) and (6.22), a specific element of the 
information matrix is 


T 1 ( qs - )T -1 ( qs ~ ) 

Mk * = Jo p i=i hi Xi(t) i R r i=i hj XjCt) i 


= h l V(k, £) h 


dt 

(6.23) 
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where V(k,£) is a qs x qs square matrix such that 

VijCk,*) = J 1 x i T Ct) T k R" 1 T £ x^t) dt (6.24) 

Matrices V(k,£) do not depend upon the input u and therefore 
they need not be computed again and again in the iterative proce- 
dure described in the next section. Notice that a considerable 
simplification results if the measurement distribution matrix, C, 
is time invariant. Then 


V ij( k,l) -|tR R' 1 T t fi.l t) xJ(t) } dt 


A Tr R- 1 T, R xx (i,j) 


(6.25) 


The quadratic constraint on the input and states may also be 
written in terms of h 

jf 1 j(T 0 x (t)) T Q(T q x (t)) + u T (t) u(t) | dt = E (6.26) 


or 


where 


H T Q h = 1 


i.. = 1 j f 1 x- 

■ij E }J 0 i 


T (t) Tj R* 1 T q x. (t) dt + AS.j 


(6.27) 


(6.28) 


6 = 1 j=j 

5 ij 0 ifij 


Case 2 

When the state and outputs are assumed multistep, the Walsh 
functions prove useful. Details of the properties of the Walsh 
functions are covered in References [56], [57], and [58]. The impor- 
tant properties for the purpose of the present application are 
summarized in Appendix B. To simplify the various derivations it 

is assumed that s = 2^, where ft is an integer. Let 4 k, 

i=0 ,1 , . . . ,s-l be the ith Walsh function in the dyadic order and 
let the input be assumed to be 
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u(t) = H $ (t) 


(6.29) 


$(t) is a vector of Walsh functions from 0 to s-1. It is shown 
in Appendix [C] that equations similar to (6.23) and (6.27) can be 
obtained for information matrix and quadratic constraint respec- 
tively . 

In the next section, we use the expressions for the informa- 
tion matrix and the quadratic constraints together with some 
theorems on certain properties of the optimal input leading to a 
computationally attractive algorithm for the selection of optimal 
multistep inputs. 


Algorithm description : Criteria for input design have been 

discussed in previous literature [59] Here, we consider two of the 
many criteria which are useful. They are: (1) A-optimality (mini- 
mize weighted trace of the dispersion matrix, i.e., minimize 

Tr(wm '*') where w is positive semidef inite) and (2) D-optimality 
(minimize the determinant of the dispersion matrix) . Extensions of 
the theorem and algorithms to other criteria is straightforward. 

In deriving the algorithms we consider only nonrandomized in- 
puts. Such an input is completely described by a vector h. An 
input h which satisfied the quadratic constraint on states and 
input is called a feasible input. It is known that the information 
matrix is always semidef inite . The following algorithms may be 
used for detection of inputs. 


Algorithm 1 

(1) Set i = 0 and select a feasible input h^ such that 
the corresponding information matrix is nonsingular. 


(2) Compute V 


m 

Z (M 

|k,&«l 


t C i ) 1 


W M 


(i) 


-1 


) kJl V(k,fc) to min Tr(WM _1 ) 


V = 


m 

Z (M 
,k,£=l 


(ir 1 


) ki v(k,£) to min |M 


-1 


(6.34) 


(3) Find the maximum eigenvalue X and the corresponding 
eigenvector h of the equation 


V h = Q h 


(6.35) 
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If the eigenvalue X is close to its minimum value (Tr(WM^ ^to 

min Tr (WM~^) ; and m to min |M’^| ), stop. Otherwise, proceed 
to (4) . 

(4) Find an input 

h tl+1) = /EE Cl) + /IE (6.36) 

by choosing a and 0 such that it is feasible and minimizes the 
corresponding cost function. 

(5) Change i to i+1 and return to (2) . 


The above algorithm converges to the global optimum. 

NOTE: What these theorems have achieved for us is an algorithm 

which converts a complex nonlinear problem into a succession of 
linear eigenvalue problems. Though the matrix V may be of rela- 
tively high order, it is positive semidefinite and we are only in- 
terested in its maximum eigenvalue — a problem much simpler than 
finding all eigenvalues of a matrix. 


Extension to nonlinear systems: For rotorcraft application, 

we have to extend the multistep input algorithm developed for 
linear systems to nonlinear systems. 


Selection of parameter identifying inputs in nonlinear systems 
is extremely complex. One major problem is that the inequality 
constraint on the quadratic function of the input and response. 
Equation (6.3) cannot be converted into an equality constraint. In 
fact, it is shown in Reference [60] that to best identify parameters 
associated with a local instability observed at high angle-of -attack 
flight of airplanes, a lower amplitude input is better than a high 
amplitude input. In addition, the theorems corresponding to 
theorems 1. and 2 on the properties of optimal inputs for linear 
systems cannot be proven for nonlinear systems. Here we describe 
an approximate method for input selection for the following system. 


x = f (x ,u, t , 0) 
y = h (x ,u , t ,0) + v 


(6.37) 


where v is white Gaussian noise with power spectral density R. 


Algorithm 2 

(1) Select an input u 0 (t) which satisfied the constraint of 
Equation (6.3) and gives a nonsingular information matrix . Let 
the corresponding state trajectory be Xg(t). 
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( 2 ) 


Linearize the system around Xg (t) , Ug (t) 
Ax = A(t,0)Ax + B(t,0)Au 
Ay = C(t,0)Ax + v 


to get 


(6.38) 


where 


Ax(t) = x ( t ) - x Q (t) 

Au(t) = U(t) - Ug(t) 

(3) Design an optimal input Au(t) 
tion (6.38) under the constraint 


(6.39) 

for the system of Equa- 


f j Ax(t) 
J Q \ 


+ x 0 ( t ) T Q Ax (t) + x Q (t) + 


(6.40) 


+ Au(t) + u Q (t) T Au(t) + u 0 (t)| dt < 1 

Au(t) may be obtained by a modification of Algorithm 1. 

(4) If Au(t) is small, stop. Otherwise, update u(t) and 

x(t) . 

(5) Return to (2) . 


Figure [6.1] is a flow chart of the rotorcraft nonlinear model 
input design algorithm. The algorithm is made up of Algorithm 1 
and Algorithm 2. First, a trajectory from a previously estimated 
nonlinear model is generated and linearized. A perturbation input 
is synthesized for this linear model based on Algorithm 1. The 
perturbation input is then added to the original nominal input and 
relinearization of the nonlinear model response to the total input 
computed. The iteration is continued until a minimum variance 
solution has been obtained. 
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CHAPTER VII 

DATA PROCESSING RESULTS 
CH-53 Results 


The CH-53A data processing results fall into two cate- 
gories: simulation data processing and flight data process- 

ing. For the simulation data processing phase, a mathematical 
model was simulated and the outputs contaminated with random 
measurement noise. These noisy outputs were then used for 
validation of various computer programs. The flight data 
processing phase uses actual flight data when executing the 
computer programs. 


CH-53A simulation data processing . - Simulation data pro- 
cess ing~ToY _ TEe _ CH r TlA - TonsTsTe3~oTi simulating the CH-53A 
and contaminating the outputs with random measurement noise, 
and processing this data through the SCIDNT, NLSCIDNT and 
DEKFIS computer programs. 


CH-53A simulation: A nine degree of freedom linear 

model for the CH-53A at 100 kts. was used for the simulation 
data base. This was a 14 state model with the body states: 
w, p, q, r, <j> and 0; and the rotor states: 3 0 >B ,8^ c> 


u, v, 

0 


lc 


S ls and 8 ls ■ 


There were 4 control inputs: B 


e Tall Rotor and 8 Main Rotor- and 14 out P uts 


Is : 
a. 


'Is 


x ’ y ’ z ’ 
The linear model 


0, 6 0 , 0 lc and B ls 
This model was simulated to the input- 


is 


p, q, r, p, q, r, <j> , 

shown in Table 7.1. 
shown in Figure 7.1, which is superior to the more common step 
or doublet input [Refs. 61 and 62]. The outputs were contam- 
inated with a zero mean white Gaussian measurement noise with 
variances: (units — ft., rad., sec.) 


R = diag [3*0.25, 3*0.0003, 3*0.000076, 2*0.000003, 3*0.000012] 


The data record was 5 seconds long and sampled at 100 samples 
per second for a total of 500 data points. 

SCIDNT runs on CH-53A simulation data: The linear identifi- 

cation program was used to identify 60 of the parameters in the 
model. The start-up values of the SCIDNT run were arbitrarily 
chosen as the simulation values _+25 percent. 

Table 7.2 shows a comparison of several parameters . 

(The quantities 3 ls /6 ls > 8 lc /8 ls and '^lc^ls 


Because of the extensive data processing done under this 
contract, only a portion of the overall results have been 
presented. An attempt has been made to present the more 
interesting and potentially useful information. 
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TABLE 7.1 

NINE DEGREE OF FREEDOM LINEAR MODEL FOR CH-S3A - 100 KTS - NO SAS 


5IA1E 

ntrivAi ivT^ — 
or siaie "v 

(fT/sir) 

(fl/stc) 

(n/scc) 

P 

(RAP/StC ) 

9 

( RAO/ SCC ) 

r 

(RAD/ SEC) 

4 

(RAD) 

e 

(RAO) 

p 0 

(RAD) 

*o 

(RAO/SEC) 

*lc 

(RAD) 

*lc 

(RAD/SEC) 

P ls 

(RAD) 

"is 

(PALI/ SEC) 

8 ls 

(RAO) 

*ls 

(RAD) 

°IR 

(RAD) 

°C 

(RAD) 

u. (ft/sec 7 ) 

-.0040 

.00042 

.0942 

-.126 

-9.97 

.865 

0 

-32.2 

25.1 

.0416 

-24.0 

-2.51 

-1.91 

-.49 

-4.92 

-4.92 

-.48 

26.5 

i. (rt/sec') 

.0156 

-.056/ 

.00216 

6.63 

-.659 

-167. 

3?.? 

0.054 

-10. 1 

-.3/6 

4.32 

-.210 

25.8 

.923 

12.0 

11.9 

6.67 

-.708 

w, (fl/sec 7 ) 

-.124 

-.00/1 

-.744 

-6.1? 

193. 

4.40 

0.89 

-1.44 

-218. 

.799 

-11.1 

3.78 

11.7 

-4.49 

124. 

17.1 

.534 

-269. 

P, (rad/sec 7 ) 

-.0011 

-.0I5R 

-.0099 

-.1/2 

.0170 

.443 

0 

0 

-.520 

-.01 

-4.56 

-.472 

17.8 

1.07 

.913 

3.33 

3.85 

-3.45 

9, (rad/sec*) 

-.0005 

.0023 

.0016 

-.0055 

-50 

.0045 

0 

0 

-.064 

-.0421 

3.84 

.30 

1.00 

.104 

-.352 

-.0795 

.0078 

.66 

r. (rad/sec 7 ) 

-.0016 

.0001 

-.0002 

.111 

.0594 

-.728 

0 

0 

-1.4 

-.0152 

-.83 

-.0532 

1.84 

.08 

-.0905 

-.543 

-5.17 

4.7 

(rail/ sec) 

0 

0 

0 

1.0 

0 

.044/ 

0 

0 

0 

0 

0 

0 

0 

0 

P 

0 

0 

0 

A. (rad/sec) 

0 

0 

0 

0 

1.0 

0 

n 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

f' 0 . (rad/ sec) 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1.0 

0 

0 

0 

0 

0 

0 

0 

0 

P Q , (rad/ sec 7 ) 

.0294 

.0005 

.432 

.495 

-3.9/ 

-3.96 

0 

0 

-168. 

-13.7 

19.0 

.999 

-16.6 

-.59 

-64.2 

9.6? 

-1.94 

231. 

r lc . (rad/sec) 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1.0 

0 

0 

0 

0 

0 

0 

f, c , (rad/sec 7 ) 

.0091 

-.0/0/ 

.0379 

-8.46 

-18.8 

.156 

0 

0 

89.5 

3.32 

-47. 

-19.4 

-161. 

-23.7 

-117. 

163. 

-.53 

30.7 

(rad/ sec) 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1.0 

0 

0 

0 

0 

Ij 5> (rad/sec) 

-.0644 

-.0506 

-.0945 

-15. 2 

12.1 

-.7? 

0 

0 

-8.31 

1.56 

154. 

22.8 

-48.7 

-17.6 

163. 

110. 

.29 

-84. 



TIME. SEC. 


FIGURE 7.1 CONTROL INPUT FOR CH-53A SIMULATION 




TABLE 7.2 - TYPICAL CH-53 PARAMETER ESTIMATES 
AFTER TWO ITERATIONS 


V s o 


V 6 o 


V 6 lc 

?1c /8 lc 

8 ls /8 1c 


V^lc 

8 lc /8 1s 

, 8 ls{ 8 ls 

8 o^1s 

8 1c^1s 

8 ls /8 ls 


8 </ 8 ls 
6 1 c /B 1 s 


SIMULATION VALUE 

IDENTIFIED VALUE 

F-RATIO 

-0.0048 

-0.00075 

0.33 

-0.0005 

-0.0017 

248.9 

-0.0567 

-0.069 

60.8 

0.00216 

0.046 

0.008 

-0.744 

-0.37 

414.5 

0.0016 

-0.0066 

150.1 

-0.0055 

0.057 

0.08 

0.111 

0.099 

68.2 

0.017 

0.47 

0.05 

-5.00 

-5.40 

64018.8 

0.0045 

-0.0349 

0.02 

-0.728 

-1.19 

781.6 

-168.0 

-138.15 

1747.5 

89.5 

42.80 

609.3 

-13.7 

-9.59 

494.3 

3.32 

5.89 

8.6 

0.78 

-2.29 

23.8 

19.0 

24.65 

375.5 

-47.0 

-51.14 

523.5 

154.0 

175.13 

168157.0 

0.30 

0.28 

2675.7 

1.00 

1.24 

97.95 

-161.0 

-172.84 

3750.6 

-48.7 

-29.25 

761.3 

-0.59 

-3.17 

431.5 

-23.7 

-19.92 

12022.3 

-17.6 

-22.22 

14004.9 

-4.92 

-5.39 

115.4 

12.0 

16.59 ’ 

1185.0 

124.0 

101.71 

19243.9 

0.913 

1.29 

968.9 

-0.352 

-0.387 

486.6 

-0.091 

-0.0075 

27.7 

-64.2 

-49.67 

8558.6 

-117.0 

-149.7 

55842.1 


L 3 3 










are rotor stability derivatives, e.g., S ls /8 ls = 38 ls /3B ls ) . 

Only two iterations on each of the 60 parameters were made. 
Correspondence is reasonable, with the exception of those 
with low F-ratios signifying poor identifyability . 

Table 7.3 shows a comparison of the starting and final 
poles of the transfer function of the vehicle response. The 
simulation values are also given. 

CH-53A flight data processing . - Flight data processing 
for the CH-53A consisted of: reading the flight data tape 

and preprocessing the data, and processing this data through 
the DEKFIS, OSR and SCIDNT computer programs. 

CH-53A flight data: The CH-53A flight data was provided 

by the Sikorsky Aircraft Division of UTC under contract by NASA. 
The maneuver processed was specifically flown at 100 kts for 
system identification purposes. The control input sequence for 
this maneuver is shown in Figure 7.2. The data record was 27.5 
seconds long and sampled at 100 samples per second. The 
following measurements were available: a x , a^, a z , p, q, r, 

V, 6, a, p, q, r, <j> , 0, » ^LAT ’ ^PED ’ 6 COLL 

PILOT PILOT PILOT PILOT 
(pilot control deflections), $ L ONG, 5 LAT’ 5 PED’ 5 COLL ( aux 

AUX AUX AUX AUX 

servo deflections in equivalent stick deflections - sum of 
pilot + SAS) , 8i, i=l,6 (blade flapping angle for all blades) 

and cos ^ (rotor azimuth). 

DEKFIS runs on CH-53A flight data: The CH-53A flight 

data was processed through the fuselage/gust estimator and 
the rotor state estimator options of DEKFIS. 

The primary purpose of the fuselage gust estimator is 
to provide data consistency between the measurements by 
correcting for bias and scale factor errors and eliminating 
random noise effects. A secondary purpose is to provide 
estimates of u, v, w and their derivatives since the meas- 
urements are V, 8, a, a x , a^ and a z . 

Figures 7.3a and 7.3b show the actual measured and esti- 
mated time histories of the normal accelerometer and the 
angle-of-attack vane. In botn cases, the estimated measure- 
ments track the actual measurements well. There is a well 
defined oscillation on the a-vane measurement (attributed 
later to boom natural frequency which was not included in 



TABLE 7.3. - COMPARISON OF IDENTIFIED CH-53A POLE LOCATIONS 



STARTING POLES 

CONVERGED POLES 

SIMULATION VALUES 


-10.1 + 1.1 7j 

-13.9 + 20 . 1 b j 

-11.46 + 21. 7j 

- 

-8.4 + 19. 9j 

-6.48 + 3. 54j 

-6.81 + 2.88J 


-4.39 + 10. 2j 

-4.68 + 10. 5j 

-6.26 + 10. 9j 


-0.246 + 1.33j 

-0.334 + 1 .25j 

-0.210 + 1 . 23 j 


-0.0123 + 0.0874J 

-0.0120 + 0.0920 

-0.0145 + O.llOj 


-5.59 

-4.75 

-4.41 


-1 .48 

-2.29 

-2.34 


-1.21 

-0.959 

-0.955 

— 

-0.219 

-0.285 

-0.167 


TIME, SEC. 


SAMPLED AT .01 SEC. 


FIGURE 7.2 CH-53A FLIGHT DATA PROCESSING - SI MANEUVER 
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FIGURE 7.3a _ COMPARISON OF STATE ESTIMATE AND CH-53A FLIGHT 
DATA (100 KTS) - NORMAL ACCELERATION 





FIGURE 7.3b- COMPARISON OF STATE ESTIMATE AND CH-53A FLIGHT 
DATA (100 KTS) - ANGLE- OF- ATTACK 
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the filter model). Figure 7.4 shows the scale factor error 
calibration on the accelerometers and gyros for this test. 

The comparison is based on the value of scale factor assumed 
prior to the filtering, indicating the change in scale fac- 
tor error concluded by the additional information provided 
by the flight data. Table 7.4 shows the bias error estimates 
obtained from this data processing. 

The rotor state estimator was used to process the six 
main rotor blade flapping measurements and the rotor azimuth 
measurement and estiamte the fixed system states (coning 

angle) , B-^ c (longitudinal flapping) and 0^ s (lateral 

flapping) . Initially, the trim portion of the maneuver 

(first second of data) was processed to obtain any biases on 

the blade flapping measurements. Coning, longitudinal and 
lateral flapping were approximated as constants in trim. 

Then, the entire maneuver was processed holding these biases 
fixed, and allowing 0 , 0^ c and 0^ s estimates to vary in 

time. The rotor state estimator also provided first and 

second derivatives of 0,0.., and 0, for subsequent use 

• 0 1C IS 1 

m OSR processing. 

An example of the estimation of flapping angles with 
the rotor state estimator option of DEKFIS is shown in 
Figure 7.5 for blade 5. 

OSR runs on CH-53A flight data: One of the key issues 

in rotorcraft modeling is the degree to which an explicit 
rotor model should be included in a stability and control 
simplified linear model. The objective is, of course, to 
work with the minimum order model which predicts stability 
or control requirements adequately. The model structure es- 
timation was applied to CH-53 flight test data to determine 
the significance of the rotor in explaining the vehicle 
response. Three generic models were formulated: 

(1) A six degree of freedom (6 dof) model with fuse- 
lage states (u, v, w, p, q, r, <p , 8). This is a quasi- 
static model in that the rotor dynamics are assumed much 
faster than the fuselage/body dynamics and are lumped in 
with the fuselage/body dynamics. 

(2) A nine degree of freedom ( 9 dof) with a first 
order rotor model with the nine fuselage states and three 
rotor tip-path-plane states (0 Q , 0^ c , 3 ls ). The first order 

rotor assumes a tip-path-plane but assumes only first order 
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FIGURE 7.4.- SCALE FACTOR ERRORS IN ACCELEROMETER AND GYRO 

MEASUREMENTS (FROM CH-53A FLIGHT TEST DATA, 100 KTS) 


TABLE 7.4.- BIAS CALIBRATION FOR CH-53A 
FLIGHT DATA (100 KTS) 


<f> 

0 

q ■ 

a x 

a y 

a z 

0.0089 

RAD 

0.0165 

RAD 

-0.0013 

RAD/SEC 

0.73 

FT/SEC 2 

-0.11 

FT/SEC 2 

-0.25 

FT/SEC 2 





dynamics for each rotor degree of freedom. (They do couple, 
however). This model corresponds to a rotor precession model. 


(3) A nine degree of freedom with a second order rotor 
model. Now, second order dynamics are assumed for each tip- 
path-plane degree of freedom. The additional states are the 
tip-path-plane rates (8 , 3^ , 3^ ). 


Each of these three rotorcraft models was identified for the 
augmented and unaugmented CH-53A. Referring to Figure 7.6, it 
is clear that closed loop parameters (with the SAS lumped in) may 
be identified by formulating the identification with J^pjlot as 

the control input, and y as the measured output. It is also 
possible to identify the unaugmented helicopter by using 
— AUX SERVO as t ^ ie contr °l inputs. As long as there is sufficient 


control excitation, the 6 


-AUX SERVO 


channel will be independent 


of the states being fed back. One might expect the identification 
of the unaugmented vehicle to provide better results since the 
SAS dynamics and nonlinearities (e.g., saturation) are inherent 
in the S_AUX SERVO i n P ut an ^ are not being lumped into the model. 

In addition, problems due to stick slop and hysteresis that occur 
between the pilot station and the auxiliary (SAS summing) servo 


are avoided, so that the 6 


-AUX SERVO 


input also provides a 


cleaner control input for identification, 


OSR identified models for each of the six cases are 
presented in Tables 7.5 through 7.10. It is possible to assess 
the quality of the results by observing the statistics in the 
tables. Any given run generates one column in the tables. The 
R2 (multiple correlation coefficient squared) at the bottom of 
any given column is indicative of the fit error, being the ratio 
of the regression sum of squares, to the observation sum of 
squares. The F-ratio at the bottom of any given column gives the 
total F-ratio for that regression. The F-to-remove number 
provides F-ratio for that regression. The F-to-remove number 
provides a measure of the relative identif iability within a given 
column. A high F-to-remove indicates a well identified parameter 
(e.g., M„ in Table 7.7). A low F-to-remove indicates poor 
b lc 

identif iability . Note that it is not advisable to compare 
F-to-remove values between column, unless the columns have a 
similar overall fit. 
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FIGURE 7.5 COMPARISON OF STATE ESTIMATE AND CH-53A FLIGHT DATA 

(100 KTS) - BLADE FLAPPING ANGLE 
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FIGURE 7.6 BLOCK DIAGRAM FOR CH53A SYSTEM 
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Results for three of these models are cross tabulated in 
Table 7.11 for the q dependent variable regressions. The 
corresponding regression fits are presented in Figure 7.7. Based 
on these, the model structure estimation program analyzed the 
relative significance of the three models with two criteria, the 

R2 multiple correlation coefficient and the F-ratio. Figure 
7.8 illustrates the relative significance of explicit rotor 
modeling. By both criteria, the increase in model accuracy is 
notable by including at least the rotor precession. Further, but 
less significant, improvement is obtained by the complete tip 
path plane model. Note, however, that in Table 7.11 the first 
order rotor model tended to enter the controls into the 
regression; whereas, one would expect most of the control effects 
to come through the rotor flapping terms. The second order rotor 
model does not enter these controls into the regression; rather, 
the tip-path-plane rates are entered giving an improved fit! On 
the other hand, inclusion of the rotor dynamics does not always 
lead to an improved fit. One can observe this in the X-force 

(u) and Y-force (v) equations. Although the r 2 multiple 
correlation coefficient increases, (which is generally expected 
since the number of terms in the regression has been increased), 
the total F-ratio decreases. This can be attributed to the 
independent regression variables fitting the noise and not the 
process. This is generally the case whenever there is a low 
signal-to-noise ratio resulting from insufficient signal 
excitation. 

Table 7.12 shows a comparison between the SAS on 
quasi-static model and the SAS off quasi-static model. It is 
interesting to note that the SAS off Mq stability derivative is 
positive (but with a very low F-ratio) and it becomes negative 
with the inclusion of SAS. (Mq may well be positive for the SAS 
off case sincere are dealing with a reduced order (quas i -static) 
model that has the rotor dynamics lumped in. It is difficult to 
be certain in this case since Mq is basically unidentifiable.) 
For the higher order models, the H SAS on effect manifests itself 

in the M and 8i c /q stability derivatives. 

H 

SCIDNT runs on Ch-53A flight data: a limited number of 

SCIDNT runs were attempted on Ch-53A flight data. Because of the 
high order of the model (14 states), the large number of 
parameters involved (>100) and the number of data points being 
processed (2750 time points for 19 channels) it was determined to 
be infeasible to run SCIDNT on the fully coupled 9 dof 
helicopter. (A problem of this magnitude, whereas possible with 
the SCIDNT software, would be extremely costly in terms of 
computer time, and would require a large number of iterations to 
insure proper cnvergence* . ) Hence, a high order identification 
model was abandoned and a simple 3rd order model (states 9, q, 
8i c ) was adopted, using only the longitudinal stick input 
portion of the maneuver. The modeling errors were compensated by 
including process noise in the model. Unfortunately, the 
converged parameter values were not reasonable, and proved to be 
very dependent. 


* The numerical efficiency of the current version of SCIDNT 
(1980) has been improved 20:1 relative to the version used for 
this study. ^4^ 
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TABLE 7.5 


CH-53A FLIGHT DATA - OPTIMAL SUBSET REGRESSION RESULTS 


• 100 kts, with SAS 

• Aerodynamic terms only 

• Si Maneuver 

• 6 dot quasistatic model 



u 

ft/sec 2 

r 

ft/sec 2 

w 

ft/sec 2 

P 

RAD/sec 2 

9 

RAD/sec 2 

r 

RAD/sec 2 

u, ft/sec 


-.1722/29.24 

-.07939/3.812 

.008592/2.413 


-.006343/20.30 

v, ft/sec 


-.08335/126.9 

-.Q3447/14.78 


.001572/32.05 

.003533/117.5 

w, ft/sec 


-.03696/3.937 

-.05283/5.826 

.01017/11.28 

-.01339/327.7 


p, RAD/sec 

-1.737/9.040 

-2.325/13.13 

1.542/6.084 

-1.075/123.0 


-.1045/13.36 

q, RAD/sec 

5.100/6.292 

3.467/3.846 

18.98/71.32 


-.6639/92.09 

.5666/56.02 

r, RAD/sec 

2.881/2.400 

-3.387/2.983 


1.854/29.01 

-.3331/17.80 

-.6182/57.32 

6 long* inches 

2.318/69.69 

.8789/14.36 

4.555/236.5 


-.1435/231.7 

.03488/11.96 

^lat’ * nc ^ es 


.3150/4.691 


.3841/233.5 

-.02354/24.06 

.04597/48.00 

6 coll, inches 


-.5512/7.278 

-7.295/780.4 

.06481/3.838 

.03535/18.91 


6p ed . Inches 


-2.597/22.15 


.2188/5.555 

-.08264/14.19 

-.3877/247.0 

F- RAT 10 

23.27 

30.83 

154.0 

54.64 

134.87 

76.62 

R 2 

.1459 

.3639 

.6655 

.4137 

.6660 

.5312 


_~J -_J 


J J 


J _J 


J 
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TABLE 7.6. - CH-53A FLIGHT DATA - OPTIMAL SUBSET REGRESSION RESULTS 


• 100 kts., with SAS 

• Aerodynamic terms only 

• SI Maneuver 

• 9 dof model - 1st order rotor 



u 

ft/sec 2 

V 

ft/sec 2 

HH 

P 

RAO/sec 2 

q 

RAO/sec 2 

r 

RAD/ sec 2 

Bo 

RAD/ sec 

"ic 
RAO/ sec 

8 ls 
RAO/ sec 

u, ft/sec 


-.1730/30.37 

-.06570/5.219 

.005569/4.989 

-.002059/6.290 

-.005380/16.62 


.001944/5.195 


v, ft/sec 

.01374/2.660 

-.08391/135.0 

-.02101/11.73 

-.001422/6.045 


.004069/173.7 



-.0009028/15.92 

w, ft/sec 


-.05269/7.931 

-.04706/3.349 


-.001211/3.737 

-.006328/37.31 


-.002832/19.78 


p, ft/sec 

-1.521/5.960 

-1.140/4.902 


.3336/31.26 

-.1356/52.34 

.1015/21.09 

-.09729/28.47 

-.1103/50.93 

-.3032/153.0 

q, ft/sec 

3.644/2.754 

3.495/4.032 

4.654/8.728 


.2000/14.50 

.2030/6.214 

-.4021/64.94 

-1.425/687.5 

-.1587/8.388 

r, ft/sec 

3.058/2.114 

-4.283/4.878 

-4.266/6.527 

.9565/41.03 - 


-.8988/111.5 



-.1351/6.725 

Bo, RAD 



-308.4/302.0 

-6.155/27.81 



1.099/6.032 

3.533/132.4 

2.544/28.24 

81c, RAD 

-12.44/2.755 


-50.16/50.67 


6.334/898.6 

-2.718/61.01 

-.6362/17.41 

-2.489/114.7 


81s, RAO 

-16.86/4.749 

28.10/20.83 

25.34/19.22 

29.61/1939. 

-1.113/25.64 

3.572/170.4 

-.9623/20.67 


-.9873/12.79 

5 long. Inches 

1.838/23.66 

.9946/18.70 



.04000/22.70 

-.3280/5.820 

-.04142/24.55 

-.1280/227.6 


Slat, Inches 




-.03842/6.392 

-.008685/3.152 


.01654/11.65 

.01268/11.17 

.04978/63.57 

Scoll, Inches 


-.6403/10.10 

-2.079/45.84 

.07251/10.24 

.01421/7.221 


.01907/5.053 


-.02514/6.898 

Sped, Inches 


-2.516/21.38 

-.9547/4.065 

.2806/41.46 


-.4212/326.9 


-.01899/2.629 


F-RATIO 

14.99 

33.35 

281.0 

422.9 

395.7 

87.72 

13.93 

88.27 

28.01 

R 2 

.1622 

.3822 

.8390 

.8758 

.8683 

.6194 

.1708 

.5953 

.2929 


KEY: 

Parameter value 

v 

xxxx/xxxx 
F to Remove^ 
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TABLE 7.7 


CH-53A FLIGHT DATA - OPTIMAL SUBSET REGRESSION RESULTS 


■ 100 kts. . with SAS 

• Aerodynamic terms only 

• SI Maneuver 

• 9 dof model - w/ 2nd order rotor 


KEY: 

-.Parameter value 
* xxxx/xxxx ' 

F to Remove* 



u 

ft/sec 2 

V 

ft/sec 2 

w 

ft/sec 2 

P 

RAO/sec 2 

q 

RAD/sec 2 

r 

RAD/sec 2 

Bo 

RAD/sec 2 

B 1c 2 
RAD/sec^ 

^2 

RAD/sec* 

u, ft/sec 


-.1746/32.95 

-.06034/5.002 

.005909/7.374 

-.001306/3.117 

-.005953/21.56 

-.01913/7.310 



v, ft/sec 


-.07649/114.6 

-.01818/10.23 

-.002701/28.68 


.004191/192.9 




w, ft/sec 


-.07854/17.21 

-.06330/7.707 


-.002838/26.48 

-.005714/31.26 

.08361/154.1 


.04008/24.41 

p, RAD/sec 

-2.798/21.27 

2.051/5.795 


-.08598/2.264 

-.1805/272.8 

.1724/26.89 

-.7654/39.09 

-.7521/14.88 

-1.729/38.01 

q, RAD/sec 

-13.16/25.43 

11.50/23.62 


-.5186/15.05 

-.2470/22.92 

.6857/58.83 

3.018/19.14 

-3.675/41:26 

6.122/51.76 

r, RAD/sec 


-3.627/3.657 

-5.294/11.07 

.8932/48.06 


-.9109/127.4 

-1.287/8.091 

-3.113/29.10 

1.266/4.809 

6 0 , RAD/sec 

6.861/8.872 


-10.78/67.29 


-.2493/31.14 

-.1671/4,685 

-.7367/2.993 



0j c > RAD/sec 

-13.23/45.56 

5.765/16.37 


-.4539/25.86 

-.2372/40.64 

.3968/37.25 

2.496/36.42 


4.248/74.02 

6j $ . RAD/sec 

-2.988/4.027 

5.188/18.00 


-1.148/165.0 



-1.369/25.20 

-2.928/65.06 

-.6940/2.548 

6 , RAO 
0 

30.63/5.021 


-288.3/292.0 


2.187/60.72 

-1.623/10.89 

-128.0/663.6 

-49.60/113.3 


B- , RAO 
lc 

-41.12/54.86 


-62.75/95.74 


5.279/1143. 

-2.030/53.41 

21.75/79.20 

-21.86/135.0 

24.52/25.00 

b ] s , rad 


39.14/25.75 

21.36/14.92 

27.79/2612. 

-1.444/93.95 

3.897/120.6 

4.291/8.274 


-46.11/234.9 

\o„q* lnCheS 


1.706/42.87 

-.4456/3.508 




-.3986/21.11 

-1.025/117.5 

.5862/22..' 

5 Laf ,ncheS 


-.5313/7.947 


.03299/6.041 


-.01542/3.885 


.1366/10.17 

.5628/69. 

A co„ •* nches 

* 

-1.015/23.77 

-2.041/46.09 




1.843/470.7 

.9106/99.79 


6 ped- Inches 


-2.392/20.58 

-1.297/8.347 

.2951/60.80 


-.4193/359.4 

-.3290/7.714 

-.6093/15.83 


F- RATIO 

19.40 

30.34 

288.6 

526.5 

497.6 

74.36 

60.50 

42.56 

31.67 

R 2 

.2003 

.4239 

.8551 

.9071 

.8924 

.6433 

.6129 

.4412 

.3701 


J __J J ___J 


J __J 


_J 


J J 


J 
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TABLE 7.8. 


CH-53A FLIGHT DATA - OPTIMAL SUBSET REGRESSION RESULTS 


• 100 kts, without SAS 

• Aerodynamic Terms Only 

• SI Maneuver 

• 6 dof quasistatic 



u 2 
ft/sec 

v 2 

ft/sec 

w 2 
ft/sec^ 

P 2 
RAD/sec 

9 2 

RAD/sec^ 

r p 
RAD/sec 

u, ft/sec 


-.1569/25.94 





v, ft/sec 

.01573/3.175 

-.09861/198.5 


-.005183/26.05 

.001326/22.99 

.003600/136.0 

w, ft/sec 

-.09037/10.77 

-.1203/30.61 

-..3201/177.7 


-.004692/31.34 

-.001968/3.465 

p, RAO/ sec 

-1.209/3.916 

-2.774/34.11 

2.696/29.37 

-.5123/51.96 

-.05173/7.385 


q, RAD/sec 

-5.861/7.084 


-7.306/15.20 

.6538/8.639 

.1569/6.085 

.3786/24.46 

r, RAO/ sec 




1.049/25.25 


-.2050/7.749 

'"’long aux, 
inches 

2.971/66.62 

1.416/27.58 

7.103/522.1 


-.2282/468.7 

.06784/27.16 

^lat aux, 
inches 


.6672/45.70 


.3803/755.8 

.02145/34.72 

.05630/149.2 

6 coll aux, 
inches 

-.6170/6.007 

-.5674/8.529 

-8.220/1435. 

.09366/11.97 

.05507/54.00 

-.04522/26.03 

^ped aux, 
inches 



-2.481/13.86 

.3927/20.93 

-.05997/6.232 

-.4620/262.7 

F- RATIO 

16.35 

47.56 

291.48 

148.3 

191.1 

82.74 

R 2 

.1530 

.3805 

.7631 

.6570 

.7387 

.5503 
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TABLE 7.9 


CH-53A FLIGHT DATA - OPTIMAL SUBSET REGRESSION RESULTS 


• 100 kts., without SAS 

• Aerodynamic terms only 

• SI Maneuver 

• 9 dof model - 1st order rotor 



u 

ft/sec 2 

V 

ft/sec 2 

w 

ft/sec 2 

P 

RAD/sec 2 

q 

RAD/sec 2 

r 

RAD/sec 2 

6o 

RAD/sec 

6 ic 

RAD/sec 

B ls 

RAD/sec 

u, ft/sec 


-.1769/34.09 



-.001906 



.001526/3.902 


v, ft/sec 


-.1106/234.2 


-.001137/3.267 

.0005303/9.281 

.004441/184.9 

-0004905/7.606 

-.0006978/12.79 

-.001370/31.89 

w, ft/ sec 

-.1134/10.33 

-.08100/13.06 

-.09912/14.06 


-.003352/37.12 

-.003908/14.18 

.003404/39.73 

.002154/7.307 


p, RAD/ sec 

-2.532/13.95 

-4.367/50.95 

2.669/24.76 

.3134/38.12 

-.1179/94.72 

.1394/32.95 

-.09999/44.22 

-.1397/62.95 

-.2638/181.7 

q, RAD/sec 

-5.832/6.324 


2.528/2.842 


-.09554/4.715 

.3622/21.72 


-.6268/189.7 


r, RAD/sec 



-4.006/7.383 

.5212/16.30 


-.2898/15.56 


.1008/6.210 

-.1513/9.245 

8 , RAD 

63.62/5.676 


-261.4/170.9 

-6.227/26.33 




3.824/46.05 

-2.712/33.60 

B, c . RAD 


63.25/29.76 

-72.89/84.08 


7.035/836.7 

-2.B13/31 .36 

-1.920/45.79 

-4.114/159.6 


B )s , RAD 

-62.06/16.33 

-45.12/12.18 

59.69/29.80 

31.33/920.7 


4.125/206.5 

-1.972/38.20 

-1.952/35.99 

-3.590/81.29 

61ong aux. 
Inches 

3.632/46.56 

4.030/53.10 



.09085/55.23 

-.04166/2.654 

-.1453/114.1 

-.2386/264.8 


5 1 a t aux, 
inches 

.8477/12.23 

1.259/37.70 

-.5632/10.81 

-.05912/12.62 

-.02750/128.0 


.03080/37.06 

.04059/63.64 

.07513/137.0 

''coll aux, 
inches 

-1.465/7.778 

-.9610/20.11 

-3.181/87.30 

.08012/11.61 


-.03909/15.89 

.05904/125.1 

.03114/8.744 

-.02269/6.264 

A ped aux, 
inches 



-3.085/31.07 

.2820/29.11 


-.4829/271.3 


.03795/5.355 

-.02905/2.077 

F- RAT 10 

14.48 

42.46 

304.47 

462.5 

513.1 

76.79 

30.53 

91.08 

38.44 

R 2 

.1764 

.4144 

.8496 

.8724 

.8835 

.5876 

.3110 

.6884 

.3624 


. _ J 


J 


J J .._J - — J J -_J - -J - — J — J 


J „_J J 
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TABLE 7.10. - CH-53A FLIGHT DATA - OPTIMAL SUBSET REGRESSION RESULTS 


• 100 kts. , without SAS 

• Aerodynamic terms only 

• SI Maneuver 

• 9 dof model - 2nd .order rotor 



u 

ft/sec 2 

V 

ft/sec 2 

w 

ft/sec 2 

P 

RAD/sec 2 

q 

RAD/sec 2 

r 

RAD/sec 2 

Bo 

RAD/sec 2 

S,c 2 

RAD/sec 

B ls 2 
RAD/sec 

u, ft/sec 


-.1817/37.35 


.004233/3.545 

-.001295/3.115 


-.01118/2.823 

.02233/7.334 


v, ft/sec 


-.09508/195.2 

-.009164/2.127 

-.003159/30.86 

.0005216/9.925 

.004771/205.4 



-.01333/30.08 

w, ft/sec 

09864/ 7. 880 

-.1606/47.91 

-.08390/10.75 


-.003514/40.37 

-.004621/12.63 

.1186/310.1 

.04956/35.28 

.02797/13.76 

p, RAD/ sec 

-5. 465/46. 24 

-1.454/6.977 

.8080/2.453 

-.1301/5.128 

-.1265/114.9 

.2042/43.32 

-.7297/24.86 

-.9488/27.40 

-2.521/107.7 

q, RAD/ sec 

-16.75/41.36 

2.997/2.441 


-.6353/19.19 

-.2358/21.15 

.7388/71.55 

5.486/139.0 

1.778/9.521 

2.826/17.52 

r, RAD/sec 



-4.029/8.039 

.4729/15.32 


-.3446/24.23 

-.7296/4.214 

-1.471/11.17 

.8181/2.494 

B Q , RAD/sec 

5.011/4.799 


-11.23/59.66 


-.2222/25.85 


-1.137/7.904 

-1.736/12.01 


Bj c> RAD/sec 

-16.70/71.39 

9.542/36.58 



-.1825/23.59 

.5920/81.97 

1.961/23.92 

-1.223/6.067 

2.442/19.41 

6 1$ . RAD/sec 

-7.351/22.44 

3.421/8.511 

-1.727/2.861 




-1.200/18.90 

-3.440/101.2 

-2.294/29.92 

B 0 , RAD 

143.8/23.87 


-241.4/160.1 


1.7323/39.09 

-2.873/10.69 

-142.7/752.9 

-48.33/56.35 

7.777/2.560 

B ]c . RAD 

-66.48/38.14 

89.94/55.01 

-87.31/57.51 

-3.429/12.86 

5.139/1115. 


9.844/9.475 

-50.99/165.7 

22.93/23.33 

P ls , RAO 

-99.64/37.03 


25.33/20.31 

25.53/678.0 


4.935/75.37 

10.07/11.97 

-14.20/15.53 

-90.04/383.2 

a, 

long aux, 
inches 


6.583/88.21 

-.9083/2.976 

-.1977/15.65 


.1073/31.97 

-1.415/81.32 

-2.911/224.5 

.9066/17.54 

*lat aux, 
inches 

1.732/42.64 

.3774/13.89 


.06482/15.10 

-.02470/108.7 

-.01776/3.699 

-.07459/2.509 

.3473/35.46 

1.129/217.1 

A .. 

coll aux, 
inches 

-1.389/8.648 

-1.885/5716 

-2.822/61.07 

.08207/19.52 


-.04923/7.590 

2.198/587.0 

1.357/145.8 


5 ped aux, 
inches 



-2.975/27.84 

.2844/34.41 


-.5070/327.7 

-.2679/3.969 

-.5812/12.18 


F-RATIO 

18.37 

40.36 

281.5 

404.9 

460.0 

76.04 

68.90 

45.98 

43.40 

R 2 

.2730 

.4521 

.8628 

.9076 

.8951 

.6295 

.6593 

.5636 

.4923 






























TABLE 7.11.- OSR RESULTS FOR THREE LINEAR MODELS 


MODEL 

DER I V AT I 

6 DOF 

9 DOF 

1st ORDER ROTOR 

9 DOF 

2nd ORDER ROTOR 

M u , l/(ft*sec) 

- 

-.00206/6.290 

-.00131/3.117 

M y , l/(ft*sec) 

.00157/32.05 

- 


M w , l/(ft*sec) 

-.01339/327.7 

-.00121/3.737 

-.00284/26.48 

M , 1/sec 

- 

-.13561/52.34 

-.18051/272.8 

M , 1/sec 

-.66386/92.09 

.20003/14.50 

-.24699/22.92 

M r , 1/sec 

-.33308/17.80 

- 

- 

V l/sec 2 

- 

- 

2.18694/60.72 

M 8 , c , l/sec 2 

- 

6.33428/898.6 

5.27916/1143. 

H 61s , l/sec 2 

- 

-1.11285/25.64 

-1.44412/93.95 

V l/sec 


- 

-.249272 

M B1c’ 1/sec 

- 

- 

-.23720/40.64 

M 0 ls* 1/sec 

- 

- 

- 

M JLONG’ l/(sec 2 *in) 

-.14347 

.03999/22.70 


M 6LAT’ V(sec 2 *in) 

-.02354/24.06 

-.008685/3.152 

- 

o 

M 6COLL’ 1/ ( sec * in ) 

.03535/18.91 

.01421/7.221 

- 

M 6PE d> l/(sec 2 *in) 

-.08264/14.19 

- 

- 

R 2 

.666037 

.868342 

.892402 

F-Ratio 

134.87 

395.73 

497.63 


• CH-53A Flight Data, 100 KTS 
t SAS On 
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MICH ACCaCKMUW, kAD/SK' 


1 


0 I.Q 2.0 3.0 4.0 S.O a.0 

TIME. SEC. 


0 1.0 2.0 3.0 4.0 S.O 6.0 

T!H€. SEC. 


0 1.0 2.0 3.0 4.0 S.O 6.0 

TIME, SEC. 


(a) 6 DOF Quasi-static Model (b) 9 DOF With 1st Order Rotor Model (c) 9 DOF With 2nd Order Rotor Model 

FIGURE 7. 7.^ COMPARISON OF KALMAN FILTERED FLIGHT DATA WITH THE REGRESSION FIT ON 
PITCH ACCELERATION FOR A CH-53A AT 100 KTS , SAS ON 


MULTIPLE CCRREUTIC.1 
COEFFICIENT CRITERION 


F-RATIO CRITERION 


FUSELAGE MOCEL J 

FUSEUGE WITH ROTOR J 
LAG MODEL 

FUSEUGE WITH TIP PATH 
PLANE MODEL “ 


LFUSEUGE WITH TIP 
PATH PUNE MOOEL 

FUSELAGE WITH ROTOR 
LAG MODEL 
FUSELAGE MOOEL 


£ FIGURE 7.8.- FLIGHT DATA IDENTIFIED CONTRIBUTIONS OF ROTOR MODEL TO OVERALL "FIT" 





TABLE 7.12.- OSR RESULTS WITH AND WITHOUT SAS 


MODEL 

DERIVATIVE"" — 

SAS ON 

SAS OFF 

M , l/(ft*sec) 

- 

- 

M v , l/(ft*sec) 

.00157/32.05 

.00133/22.99 

M w , l/(ft*sec) 

-.0134/327.7 

-.00469/31.34 

M p , 1/sec 

- 

-.0517/7.39 

M , 1/sec 

-.664/92.09 

.1569/6.09 

M r , 1/sec 

-.3331/17.80 

- 

MgO’ lyrsec 

- 

- 

"eic- 1/secZ 

- 

- 

M 3 l s ’ 1/sec 

- 

- 

V ,/SeC 

- 

- 

M B1c’ 1/sec 

- 

- 

M Bls’ ,/sec 

- 

- 

M 6L0MG’ V (sec^in) 

-.143/231.7 

-.2282/468.7 

M 6LAT s 1 /(sec 2 *in) 

-.0235/24.06 

.0215/34.72 

M 6C0LL’ 1/ ^ sec * in ) 

.0354/18.91 

.0551/54.00 

M 5PE d» 1/(sec 2 *in) 

-.0826/14.9 



-.05997.6.23 

R 2 

.666037 

.73865 

F-Ratio 

134.87 

191.13 


• CH-53A Flight Data, 100 KTS 

• 6 DOF Linear Model 


5 


PARAMETER VALUE 
XXXXX/XXXXX 
F TO REMOVE 


5 
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on the level of process noise chosen. Further attempts to iden- 
tify the process noise power spectral density were inconclusive. 
The moral of this effort is that it is dangerous to rely on low 
order models when significant higher order dynamics are present. 


Bell 609 Rotorcraft Results 

Parameter identification results were generated with linear 
SCIDNT for the Bell 609 rotorcraft, using actual in-flight data 
as input. This flight data was acquired from Bell Helicopter 
(Textron) on a magnetic tape which contained control and 
measurement time history redords of longitudinal maneuvers 
performed on May 23, 1973. The particular maneuver data used in 
this identification example was #4158, which consisted of a 1 
second aft pulse on the stick (longitudinal cyclic input J. This 
particular maneuver was chosen because there was very little 
lateral stick motion. 


A fourth-order linear model was used in the parameter 
identification program. This model is 


x(t) = Fx(t) + Gu (t) 

y(t) = Hx (t) + v ( t ) 

E[v(t)v T (r)] = R6(t-r) 

Assigning the states u(t), w(t) , q(t), and 0(t), long- 
itudinal velocity, vertical velocity, pitch rate, and pitch 
attitude to x and y, and assigning the controls B lc (t) 

and 1. to u, the model can be written as follows: 

-g cos8 0 

-g 5in8 0 
0 

0. 



u(t) 


w(t) 


q(t) 


9(t) 


d 

at 


u(t) 

w(t) 

q(t) 

9(t) 

' X 3 : 

Z B : 


u 

0. 


M 

w 

0. 


X -w 
q 0 

Z *U 
q o 


q 

l. 


"b, (ty 
1. 


0 0 
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u (t) 
m 


i — 1 

0 . 

0 . 



0 . 

1 . 

0 . 



0 . 

0 . 

1 . 

9 m^ 


0 . 

0 . 

0 . 


0 . 


u(t) 


n u Ct) 

0 . 


w(t) 

+ 


0 . 


q(t) 


n qU) 

1 . 


9 (t) 


n 0 (t) 


The time history measurements used were: V , total 

velocity, a , angle-of-attack, q m , pitch rate, and 

0 , pitch attitude. V and a were resolved into u, 
m r mm m 

and w , and after proper units conversion, the tape was 

read by SCIDNT. A total of 287 points, spaced at a .125 
second sample interval, were used (every 64th point was used). 
The nominal airspeed of the maneuver was 60 kts (100 ft/sec). 

The SCIDNT run was started with most of the unknown 
parameters in F and G as zero. The F and G matrices 
were initially input as follows (with parameters to be iden- 
tified starred). 


‘ * 

* 

* 


0 . 

0 . 

0 . 

-32.17 

* 

* 

* 


0 . 

0 . 

100. 

. 337 

* 

* 

* 


0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

1 . 

0 . 

* 




0 . 

0 . 



* 

* 



0 . 

0 . 



* 

* 



0 . 

0 . 



0 . 

o._ 




152 



After 27 iterations, linear SCIDNT calculated the follow- 
ing F and G matrices : 

Final F = 


- 6* 

59x10 0 

-1.83xl0' 6 * 

-31.7* 

-32.17 

157* 

4.93xl0' 7 * 

66.9* 

.337 

00567* 

_ q * 

-9.0x10 y 

-3.83* 

0 

0 

0 

1 . 

0 


Final G = 

.116* - 5.41* 

.837* -40.3* 

-.0244* 1.16* 

0 . 0 . 


The final F matrix had the following eigenvalues : 
-3.8, 0.0, and -1.73 +_ j0.22. 


The particular estimates from these matrixes, as well as 
the standard deviations and F-ratios are presented in Figure 7.9. 

It can be seen from the F-statistics that the pitching 
moment coefficients should be regarded as the most accurate para- 
meter estimates. This is to be expected, as longitudinal stick 
motion primarily excites pitch dynamics. Other parameters, such 
as X u , X w , Z w and M w are basically unidentifiable as is 

indicated by the low F-ratios. 


Since the initial value of the control was not subtracted 

out of the data, and bias terms X . Z . and M_ are estimated, 
’ 0 0 0 

an approximate value of the control bias can be calculated as 


B 


lc. 


= -<x + z; 


m^/OS 


lc 


Z R Z o + 
B lc ° 


M 


B, M o> 
lc 


48.21 
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The estimate of the control bias is fairly good, that is, 48 
percent, compared to an actual value of 50 percent. 

In conclusion, careful preliminary analysis of the flight 
data and the modeling should be made prior to any SCIDNT run. In 
particular, it is recommended that DEKFIS be used to verify 
flight data consistency and OSR be used to assure proper 
modeling. In this case, it was attempted to bypass these steps 
with obvious consequences. 


RSRA Results 

Results for the rotor systems research aircraft (RSRA) 
consist of simulation data processed through the optimal subset 
regression (OSR) computer program 

The simulation data was provided by NASA Langley Research 
Center and was generated by the Nonlinear RSRA simulation 
mathematical model [Ref. 63]. This model was simulated to the 
control inputs shown in Figure 7.10. These inputs are superior 
to the more common step and doublet inputs usually used for 
identification [Refs. 61 and 62]. Unfortunately, the magnitude 
of these control inputs was exceedingly large, and the unstable 
simulation model diverged over the length of the 25.5 second 
maneuver. Nevertheless, the data were processed through OSR to 
show its flexibility in breaking out the rotor force and moment 
contributions. Whereas the results alone are of questionable 
value, the techniques investigated should prove to be quite 
useful when processing actual RSRA flight data. 

Three different models were identified from this data. 

These were: 

(1) 12 dof with 2nd order rotor; (2) rotor hub force/moment 

breadkdown; (3) 6 dof with rotor hub force/moment contributions. 

The first model is a 12 degree of freedom (dof) model with 6 
body dof's, 3 rotor flapping dof's and 3 motor lagging dof's. 

Each rotor flapping and lagging degree of freedom is of second 
order. The resulting OSR identified model is shown in Table 
7.13. Because of the large nonzero mean excursions in the 
independent variables (i.e. ,u,w,q) , the parameters identified in 
the linear model are not physically realistic. An attempt was 
made to identify some nonlinear (polynomial type) parameters. 

This is shown in Table 7.14. 

The second model is a rotor hub force and moment breakdown 
model where the hub forces and moments were treated as dependent 
variables and the 12 dof rotor (only) model was identified. 

Table 7.15 shows this breakdown. 
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The third model identified is a 6 dof body/f useiage model 
where the rotor hub forces and moments are treated as control 
inputs (i.e., independent variables). The results are shown in 
Table 7.16 

Note that in many cases quite reasonable regression fits and 
F-ratios were achieved. These results can be deceiving, 
however. Since the sum of squares of the signal was large, the 
residual error will be proportionately large for a given multiple 
correlation coefficient (R2). These residual errors are 




Standard 


Parameter 

Estimated Value 

Deviation 

F-Value 

X u 

3.59x10'® 

1.62xl0' 3 

4.93x10'®* 

Z u 

-.157 

.0482 

10.5 

H u 

.0057 

.0005 

126. 

X w 

-1.83x10'® 

7.78xl0~ 3 

5. 56x1 0' 8 * 

Z w 

4.93xlO' 7 

6.73xlO' 3 

5.36xl0' 9 * 

M w 

-9.0xl0' 9 

1 .48xl0' 4 

3. 69x1 0' 9 * 

x 9 

-31.65 

1.84 

29.4 

Z q 

-33.08 

33.06 

4.09 

M q 

-3.83 

.336 

130.0 

V 

.115 

.023 

25.2 

z °,« 

.837 

.160 

27.26 

m r 

B lc 

-.0244 

.00171 

20.3 

X o 

-5.42 

1.14 

22.6 

Z o 

-40.4 

7.63 

27.9 

H o 

1.16 

.0818 

201.1 


* Not Identifiable 


FIGURE 7.9- PARAMETER ESTIMATES AND STANDARD DEVIATIONS 
FOR BELL MODEL 609 AT 60 KTS . 
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n 


not negligible, and substantially more nonlinear modeling is 
required to yield reasonable numbers. On the other hand, if 
too many additional independent variables are added, the number 
of independent variables can equal or exceed the number of data 
points in the regression. For the supplied sampling period 
(T=0.2 sec), it is neither reasonable to substantially increase 
the number of independent variables to model the nonlinearities 
nor to reduce the length of the data segment to maintain the 
linearity of the data. 


UH-1H Results *”} 

| 

The UH-1H results consist of control input design and 
flight data processing through the SCIDNT computer program. ^ 

Initially, a series of stepwise control inputs were designed 
for implementation in the UH-1H flight test program based on 1 

a UH-1H linear model. These inputs were subsequently dupli- 
cated by the pilot during the flight testing. The conven- “1 

tional doublet inputs were also used in the flight tests. / \ 

The final step was to process the data through the SCIDNT 
parameter identification program and get improved linear — , 

models. j 

Input design for UH-1H flight testing . - To illustrate 
the sequence of computations and type of inputs, the applica- 
tion of the input design algorithm to a low order simulated 1 

linear model of a UH-1H helicopter (no stabilizer bar) at 60 kts 
was computed. Only the longitudinal dynamics were considered, — 

and thus the quasi-static rotor assumption is evoked. The 
assumed model, based on C-81 [Ref. 64] computations, is shown in 
Figure 7.11, where five measurements are assumed (u, w, q, 0, 
a z ) and with measurement random noise power spectral density j 

(units ft, rad, sec) 

R = diag [4. ,4. , 10" 4 , 10 -4 , 4 . ] ~j 

» 

The two controls, collective stick S^OLL and longitudinal 

cyclic stick ^lqng are assume ^ applied individually. It is j 


desired to design inputs which maximize the accuracy of the 
parameters underlined in the matrices (10 parameters) . For 
this input, 8 seconds of total input time are specified, with 


only 1 

second steps 

allowed. For both 

inputs, ^coLL anc ^ 

1 

6 L0NG’ 

mean square 

control excursions 

are limited: 

1 


f 8 2 

J 100 6 L 

0 

dt = 1.0. 


1 

1 
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1 1 

15 20 



TIME, SEC 


FIGURE 7.10.- RSRA SIMULATION CONTROL INPUTS 
(STICK DEFLECTION IS IN PERCENT 
OF TOTAL TRAVEL RELATIVE TO THE 
TRIM STICK POSITION) 
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TABLE 7.13. - RSRA SIMULATION DATA - OSR RUNS - LINEAR - 12 DOF WITH 

2ND ORDER ROTOR 


„t' 7 ’•* — F to Remove 

Coefficient Value 



FT/SEC 2 

FT/ SEC 2 

FT/SEC 2 

P , 
RAD/SEC 2 

RAD/SEC 2 

r 

RAD/SEC 2 

u, ft/sec 

-.108938/659./ 

- 

-.0982755/55.82 

- 

- 

-.0001866/42.16 

v, ft/sec 

- 

.780321/582.7 

.527098/36.42 

.0040090/11.15 

- 

.0025545/61.72 

w, ft/sec 

- 

-1.32923/134.1 

- 

- 

- 

-.00131511/36.62 

p, rad/sec 

.792670/17.20 

- 

- 

- 

- 

-.0424003/45.91 

q, rad/ sec 

- 

-10.7371/5.915 

136.726/266.3 

- 

- 

- 

r, rad/sec 

-38.5933/456.2 

-154.652/403.9 

78.9239/17.53 

- 

-.500467/161.4 

-.309070/23.53 

4>, radians 

-2.75896/660.8 

- 

- 

- 

.0292561/119.6 

- 

4, radians 

-35.8445/1623. 

- 

- 

- 

-.0613266/15.35 

- 

fi, rad/sec 

- 

-38.0693/23.41 

- 

2.33184/35.98 

.253363/14.37 

-.154603/6.503 

B lc , rad/ sec 

-6.37844/75.98 

44.5946/23.29 

16.9444/19.89 

-1.64247/180.6 

-.15317/45.97 

.137973/38.46 

6 )s , rad/sec 

27.5079/72.46 

-207.816/376.6 

-190.027/133.4 

10.2840/439.9 

1.66507/265.4 

-.351517/20.73 

V rad 

- 

495.536/69.57 

- 

- 

-2.02214/170.9 

- 

g )c . rad 

- 

135.652/149.2 

125.448/23.94 

- 

- 

- 

fi| s . rad 

- 

- 

- 

- 

- 

.430655/181.5 

C , rad/ sec 

- 

- 

- 

- 

- 

- 

C )c . rad/sec 

- 

- 

- 

- 

- 

- 

f, ls> rad/sec 

- 

- 

- 

- 

- 

- 

V rad 

25.5739/98.12 

- 

-149.185/69.88 

- 

.631489/52.49 

-.662435/105.6 

C, c . rad 

- 

- 

- 

- 

' - 

- 

V rad 

.0849889/30.03 

- 

- 

- 

- 

.0014917/17.53 

fi LONG’ inches 

- 

- 

- 

- 

- 

- 

d LAT’ ’ nc ^ es 

- 

-1.66441/85.07 

-.569339/33.32 

- 

.0092254/232.6 

- 

» Inches 
inches 



' 


" 


F-Ratio 

4631.7 

179.25 

46.271 

181.12 

96.694 

82.9 

n 7 

.997177 

.938729 

.798175 

TS548T3 

.891961 

.687149 


J J J ] J _ J J J 1 . . J J ] \ .... J 
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TABLE 7.13. -CONCLUDED 



RAD/ SEC 2 

®lc 

RAD/SEC 2 

®ls 

RAD/SEC 2 

*o 

RAD/SEC 2 

e ic 

RAD/ SEC 2 

'is 

RAD/ SEC 2 

u, ft/sec 

- 

- 

- 

-.00135146/93.50 

- 

.0005515/83.09 

v, ft/sec 

- 

- 

- 

- 

- 

- 

w. ft/sec 

- 

- 

- 

- 

- 

- 

p, rad/ sec 

- 

- 

- 

- 

- 

- 

q, rad/sec 

- 

- 

- 

.487964/12.41 

- 

-.161654/19.65 

r, rad/sec 

- 

- 

- 

- 

-.111215/7.750 

- 

<h rad 

- 

- 

- 

- 

- 

- 

6, rad 

- 

- 

- 

- 

- 

- 

B 0 , rad/sec 

1.24137/14.70 

-3.70532/50.43 

-1.1000/61.52 

- 

- 

- 

e ]c . rad/ sec 

.615372/35.19 

-2.49783/255.3 

•-.506508/120.4 

- 

-.479954/1009. 

.215715/163.4 

6, s> rad/sec 

1.66034/21.50 

-3.05377/27.76 

- 

- 

- 

- 

V rad 

- 

- 

-.532184/28.29 

- 

- 

.658497/46.39 

B, c rad 

- 

- 

-1.1689/56.93 

- 

- 

- 

P 1s . rad 
C 0 . rad/sec 


_ 

-.458404/33.76 

-3.89073/20.94 

.331959/53.39 

“ 

C, c . rad/ sec 

3.36914/11.45 

- 

- 

- 

3.03866/134.7 

- 

t, . rad/ sec 

- 

- 

- 

- 

- 

- 

V rad 

- 

- 

- 

-3.92671/125.2 

- 

- 

t, c . rad 

- 

- 

- 

- 

- 

- 

'-Is* rad 

- 

- 


* 

- 

-10.5956/145.4 

6 L0NG- i,,c,ies 

- 

-.022818/44.47 

- 

- 

- 

- 

d LAT’ * nc ^ es 

-.0165573/60.87 

- 

- 

- 

-.002584/19.59 

- 

d COlL’ 1nC,,eS 

- 

- 

- 

- 

- 

- 

6 pED , inches 

- 


- 


- 


F-Ratio 

84.373 

255.34 

87.907 

47.629 


84.357 

R 2 

.775680 

.8925 

.78289 

.607678 

IlflHiiiiS 

.775646 


j j * ' / I — F to Remove 

Coefficient Value 















TABLE 7.14. - RSRA SIMULATION DATA - OSR RUNS 


- Longitudinal 

- Transgenerated u 2 , w 2 , q 2 , 

- 12 dof w/ 2nd order rotor 


V/**- 

Coefficient Value 


"F to Remove 


uw, u, q, wq terms 



u 

w 

q 


ft/ sec 2 

ft/sec 2 

ft/ sec 2 

u, ft/sec 

.0462094/45.23 



v, ft/sec 

-.150251/328.7 



w, ft/sec 



-.00222258/7.712 

p, rad/sec 


-51.4330/169.9 


q, rad/sec 

-12.0530/54.60 

50.3129/38.40 

.287329/45.96 

r, rad/sec 


90.4102/72.66 

-.479089/103.5 

u 2 , ft 2 /sec 2 

-.000282046/709.7 



w 2 , ft 2 /sec 2 




q 2 , rad 2 /sec 2 


-331.375/56.42 

-1.90977/72.51 

u*w, ft 2 sec 2 


.00720504/153.0 

-.000014654/20.38 

u-q. ft-rad/sec 2 

-.172689/132.4 

.976432/139.4 


wq. ft-rad/sec 2 

-1.16709/213.7 


.0187846/43.95 

q, radians 

-.733555/43.35 


.0260525/153.6 

<p, radians 

-27.2399/1544. 



e 0> rad/sec 




Sic » rad/sec 



-.220495/112.3 

6 1 2 » rad/sec 



1.30155/398.3 

Sq, radians 


199.449/39.26 

-1.48749/12.30 

6 ic, radians 

8.79354/47.60 


-1.01747/119.4 

B -j 5 , radians 


840.703/189.9 


; 0 , rad/sec 




; lc , rad/sec 




C1S> rad/sec 




C 0 , radians 


-153.125/74.72 

.591002/97.90 

5 l radians 




CIS, radians 




5 L 0NG> inches 


-.595814/49.57 


5lat » inches 



.00750275/42.01 

°C0LL ’ inches 



.00750275/42.01 

Cp£D> inches 




F-Ratio 

10210. 

49.731 

129.75 

_R 2 

.998718 

.809541 

.936692 
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TABLE 7.15. - RSRA SIMULATION DAT 
2ND ORDER ROTOR - P 


Coeffi- 

cient 

Value 


u, ft/sec 

v, ft/sec 

w, ft/sec 

p, rad/sec 

q, rad/sec 

r, rad/sec 
rad 

B, rad 

6 L0NG- inches 
dlftT. inches 
6 C0L l. inches 
<Sp£Q, inches 
U 0 * rad 
die. rad 
dl 5 > rad 
0 o . rad/sec 
die. rad/sec 
^ 1 $, rad/sec 
C 0 . rad 
Ci C . rad 

C, s . rad 
C 0 . rad/sec 
ClC> rad/sec 
Cis. rad/sec 

F-Ratio 

R 2 


X R 

ft/sec 2 

y r 

ft/sec 2 

-.00652639/45.39 

.02646/51.93 

-2.1530/62.34 

-.192164/79.05 

2.39862/155.0 

.030524/66.35 

1.84322/615.0 

1.28644/3.263 


-.0240146/24.90 

-1.39531/131.1 

-1.95725/122.8 

7.19282/92.32 


-2.51541/56.39 

2089.3 

383.57 



601 .07 


6 


1094.0 

.901900 


924.32 

.904162 
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TABLE 7.16. - RSRA SIMULATION DATA - OSR RUNS - LINEAR - 6 DOF WITH 

ROTOR FORCE/MOMENT CONTRIBUTIONS 


Coem- 
clent Val 

r ' I 

u 

ft/sec 2 

V 

ft/sec 2 

ft/sec 2 

P 

rad/sec 2 

9 

rad/sec 2 

r 

rad/sec 2 

K to 

lue Remove 


u, ft/sec 

-.106146/360.9 





-.000151235/2.313 


v, ft/sec 






.00917058/507.2 


w, ft/sec 



.627368/50.99 


-.00851137/341.4 

-.00211978/66.10 


p, rad/ sec 


-33.2212/144.4 



-.0302990/10.27 

.110442/52.49 


q, rad/sec 

14.3666/48.48 


110.979/132.9 





rad, rad/sec 

-22.1317/112.2 

-49.6356/27.72 

96.4448/90.56 



-1.14855/209.1 


<|>, rad 

-2.68670/286.5 

-3.38161/42.00 

4.87067/115.7 

.0320097/231.5 




0, rad 

-41.5172/1474. 

13.9795/73.52 




.0991914/16.96 


•Tong- ' nches 





-.00149956/8.895 



5 LAT’ 1nches 








a C0LL" lnches 


.493362/44.62 

576395/23.85 



-.00557196/176.6 


6p ED , inches 








X H . ft/sec 2 

3.51784/74.48 




.0743372/30.05 

.0102781/6.407 


Vr, ft/sec 2 


15.1571/68.35 




-.0367486/17.04 


Zr, ft/sec 2 


.357942/28.54 

1.40466/105.1 

| 

-.00617716/40.11 



Lr, rad/sec 2 


-40.0115/259.8 


1.10766/12460. 

.0958124/145.0 

-.0538265/32.54 


Hr, rad/ sec 2 

12.3210/12.26 


32.6719/17.82 


1.2/048/38.70 



Hr, rad/sec 2 




.522857/84.28 

-.572691/74.12 

.994565/613.1 


F-Ratio 

3840.6 

111.87 

33.763 

4596.4 

94.979 

139.06 


R 2 

.995556 

.882638 

.663244 

.991088 

.864593 

.929513 


J __ J . J 


J J . J 


J J 


-J J j 


J 
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FIGURE 7.11. -THREE DEGREE OF FREEDOM LINEAR STATE 
AND MEASUREMENT EQUATIONS FOR THE 
UH-1H HELICOPTER 

• No Stabilizer Bar 

• 60 Kts 



Figures 7.12 and 7.13 show the iteration outputs as the 
algorithm minimizes the sum of the estimated variances of the 
parameter error covariance matrix (J) . In both cases, a con- 
ventional pulse input is used to start the iteration. The 
cyclic input is seen to reduce the cost J by 60 percent. 

The collective input reduces J by 68 percent. 

A considerably more complex input design approach was 
taken for flight test. The linear math model used in this 
design was derived from C-81 simulation. It is an eight 
degree-of-freedom, 10-state model with four inputs. It is 
assumed that we have measurements of 10 independent quanti- 
ties. The states, controls and measurements are as follows: 

States (x) : u, w, q, 0, v, p, r, <p , a^, b-^ 

In P - ut - .. ? . — ( hJ. • 5 L0NG’ 5 LAT’ 6 PED’ 5 coll 

Measurements (y) : u, w, q, 8, v, p, r, <f> , a^, a z 

The state and measurement equations are 
x = Fx + Gu 

y(k) = Hx(k) + Du(k) + v(k) 

Nominal values of F, G, and H matrices used in input de- 
sign are shown in Table 7.17. Note that we have continuous 
state equations with discrete measurements. v(k) is a noise 
term which is assumed zero mean white Gaussian with covariance 
(units ft., sec., rad). 

R = diag [4.0, 4.0, 0.0001, 0.0004, 4.0, 0.0001, 
0.0001, 0.0004, 0.0004, 4.0] 


There are about 80 parameters of interest. Since a single 
input could not possibly excite the system to identify all the 
parameters simultaneously, 10 individual sequences were de- 
signed. Each input sequence was designed to maximize estima- 
tion accuracy of only 10 parameters. Some important parameters 
are included in two or more sequences. Table 7.18 lists the 
parameters for which each of the 10 inputs is designed. The 
inputs are shown in Figure 7.14. Each input consists of a 
sequence of eight steps each of 1 sec. duration and has been 
designed to minimize total input energy (which in turn also 
restricts output energy) . 




FIGURE 7.12 UH-1H INPUT DESIGN - LONGITUDINAL STICK 

CONTROL, 6 L0NG 


[TERATlOfi 

8 STEP INPUT 

J 

LARGEST 

EIGENVALUE 

Y 

1 

5 COU 


L_ » t 

68051. 

104385. 

.77643 

2 

A CnLL 



46599. 

47168.3 

.99022 

3 

5 coll 

^ - 

46209. 

46229.6 

.99973 

4 

4 coll 



46195. 

46196.5 

.99998 

5 

4 COLL 

^ — ► t 

46194. 




FIGURE 7.13 UH-1H INPUT DESIGN - COLLECTIVE STICK 

CONTROL 6 CO ll 


165 




























166 



1 I I I I I 1 1 1 1 1 11 1 1 ) 1 1 1 

Table 7.17 

Eight Degree of Freedom Linear State And Measurement 
For the UH-1H Helicopter 

• No Stabilizer Bar 

® 60 Kts 



F 





1 2 


1 

-l,390oE-C2 6.3300L-62 

-.7 . 

c 

9. 9000t-09 -7. 5E60E-01 

l. 

3 

2.2150E-C3 -3. 9700E-03 

-3. 

9 

6. 0. 

l.i 

5 

.l*?2.90|-0i_^2 t ?7?ot-rC2 


6 

1 .95306-03 -3. 17PCE-03 

1. 

7 

-2,35t»iE-L3 -4, 6670E-02 

1* 

6 

C. 0. 

0. 

9 

2. 9800E-03 7. 8SCOG-03 

-l.i 

10 

-9.1200E-03 5.0300E-03 

3* 


o. 

0. 

Of. 

0. 

0*. 

a. 

Q* 

a. 


_ 1.2800E-02 
-2.39C0E-D2 
-9*10001-09 

6* 

-s.3eooe-02 

9.919QE-02 

o. 

1.7900E-03 

-5.9960E-03 


5,2150fc-0l 

-i.389UE*0b 

-7,01006-02 

0. 

6,64?C|;*90 

-7.9630E-01 

3,36006-01 

l.G0J0t*U3 

1.2360E-02 

-i.C3006*0<> 


-1.8530E-01 

1.85316*93 

-1.92906-92 

0. 

I. 7059E*00 
-1,3770E*09 

7.21 33E-92 

J. 0400E-02 
-1.1550E-01 


§ 

0, 

0. 

0. 

0. 

9.2b?DEt01 

8. 

0. 

8. 

P. 

0. 


9 

-3,18006*01 

1.01876*01 

5.18706*00 

6 . 

1.15.576*01 

1.35736*01 

-6.27906*00 

0. 

-l,3955t*01 

6.56706*00 


10 

8.97006*00 

2.13106*00 

-1.96336*00 

0. 

2.07176*01 
x. 53016*01 
-2.73906-01 
0. 

-6.66506*00 

-1.60366*01 


G 


12 3 9 

1 1.2630L*CC -1.967GE-01 -3.266CE-01 -9.08C0E-02 

2 -1.323 U £+0i 3, 5630E+C0 -2.9700E-02 -9.38661-02 

. 3 -1.106OC-01 -7. 9C00E -C 3 i. 3310^-02 -1.890 0E-02 

9 c. b. c. 6. 

5 -3.221o£-pi 9.36906-01 2.15006-01 1,67816*00 

6 -3.99061-1)1 2. 93e0L-C l 2.78616-01 2.68706*06 

7 l?9l30|-tl_ 3. 1300E-C2 -1.96606-01 -1,60101*00 

8 0. 0. 0. 6, 

9 _ 1 .25C0E-01 -9. 993pt-C 1 ItBS^E-Ol 2,2900E-02 

10 ' — 6. 26JC0£— |i 2 2.3i70l-0i 3.70906-61-1.15001-02 



1 

2 

3 

9 

5 

6 

7 

8 

9 

10 

1 

i.o6ooe*6o 

do 

0. 

0. 

0. 

0. 

0. 

0. 

6. 

0. 

2 

0. 

1.C6006 ♦( 0 

0. 

0. 

C it... 

0. 

0. 

o* 

0. 

0. 

3 

0. 

0. 

i.ouoct*oo 

0. 

0. 

0. 

0. 

0. 

6. 
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TABLE 7.18. - INPUTS AND THE CORRESPONDING PARAMETERS FOR WHICH 

THEY ARE OPTIMIZED 


INPUT #1 

INPUT #2 

INPUT #3 

INPUT #4 

INPUT #5 

INPUT 16 

INPUT #Z 

INPUT #8 

INPUT #9 

INPUT #10 

6 LONG 1 

''LONG 2 

6 L0NG 3 

\at 1 

{ LAT 2 

6 LAT 3 

6 PEDAL 1 

6 PEDAL 2 

6 COLL 1 

6 C0LL 2 

H 

Y 

X 

L 

M 

X 

Y 

x„ 

Z„ 

Y w 

U 

U 

U 

V 

P 

V 

V 

V 

W 

W 

M 

L 

z 

L 

N 

z 

N 

z 

m. 

L w 

u 

w 

u 

V 

P 

V 

V 

V 

w 

w 

H 

L 

L 

L 

M 

Y 

L„ 


M 

N, 

w 

q 

u 

P 

P 

V 

V 

V 

w 

w 

H 

L 

N 

L 

M 

N 

L 

M, 

z . 

L, 

q 

q 

u 

P 

r 

V 

r 

V 

al 

w 

M 

Y . 

L 

L 

L , 

M 

N 

H 

M , 

L 

q 

al 

u 

r 

al 

V 

r 

r 

al 

q 

X al 

L al 

"q 

Y bi 

x bi 

L v 

L r 

M r 

M al 

L al 

Z al 

N al 

X al 

L bi 

z bi 

Y bl 

Y 6PED 

X 6PED 

X 6C0LL 

Y 6C0LL 

M al 

L al 

X 6LONG 

N bi 

M bi 

Y 6LAT 

L 6PE0 

Z 6PED 

Z 6COLL 

L 6C0LL 

M a1 

M bl 

Z 6L0NG 

L bi 

M bi 

L 6LAT 

N 6PED 

M 6PED 

H 6COI.L 

N 6C0LL 

M 6L0NG 

L 6LONG 

M 6LONG 

L 6LAT 

M 6LAT 

N 6LAT 

L 6PE0 

M 6PED 

M 6C0LL 

L 6C0LL 


J 


. J 1 J 1 1 1 . ] 


] . J ;_*j 




UH-1H Flight Data Processing .— Flight data processing for 
the UH-1H consisted of: reading the flight data tape, editing 

and preprocessing the data, and processing the data through the 
SCIDNT computer program. (Due to the limit imposed on computer 
resources, it was not feasible to run DEKFIS and OSR on this 
data . ) 

UH-1H flight data: Flight data for the UH-1H V/STOLAND 

helicopter was provided by NASA Ames Research Center. This data 
was from Flight 68 and dated February, 1978. The maneuvers 
processed were specifically flown at 60 kts for system identifi- 
cation purposes using the optimally designed inputs discussed 
previously. Conventional doublet inputs were also flown for 
comparison. The pilot's rendition of the optimal longitudinal 
stick input ^LONGj-* s ^ own i- n Figure 7.15. The optimal 

input sequence is repeated three times in succession. Figure 
7.15 also shows the collective stick motion during this period. 

The optimal collective stick input analyzed is likewise presented 
in Figure 7.16. Optimal collective stick input 6 is 

COLL2 

repeated three times at the beginning of the maneuver and once 
towards the end. Longitudinal and collective doublet inputs 
are shown in Figures 7.17 and 7.18 respectively. The four 
data records processed (each corresponding to the inputs of 
Figures 7.15 to 7.18 were 50 seconds long with a sample period 
of 0.0504 second (for a total of 990 data points). “The following 
measurements were utilized: ^XOTAL* 0 > a x , a z , and 

6 coll. 

SCIDNT runs on UH-1H data: The UH-1H flight data was pro- 

cessed using a linear 3 dof, 4th order longitudinal model. (A 
simple longitudinal model was chosen in order to conserve com- 
puter time.) Since all the maneuvers were flown at the same 
flight condition, it was decided to process a longitudinal stick 
input and a collective stick simultaneously to obtain a single 
identified linear model. Three separate longitudinal-collective 
stick maneuver combinations were considered. These cases are 
shown in Table 7.19. These cases produced three separate final 
estimated parameter models: a model identified using only optimal 

inputs, a model identified using only doublet inputs, and a model 
using the combination shown for Case 3 (to be discussed further 
below) . 

The initial parameter estimates were based on C-81 computa- 
tions for a 3 dof longitudinal model without stabilizer bar stabil 
ity augmentation. The F,G,H and D matrixes for this model are 
shown in Figure 7.19, with states, controls and outputs as follows 

States : u, w, q, 8 

Control Inpu t s : 6 L0NG) S C0LL 

, h, q, 8, a x , a z 


Outputs : V TOTAL 
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This initial parameter model is shown simulated to the pilot's 
optimal longitudinal stick inputs (of Figure 7.15 in Figures 
7.20(a) - (f) . This simulation ("B" trace) is shown compared to 
the actual UH-1H response time history ("A" trace). Large 
discrepancies can be noted. The same model is likewise simulated 
to the pilot's optimal collective stick input (of Figure 7.16 
in Figures 7.21(a) - (f) . 


For Case 1, a new set of parameters was obtained by running 
SCIDNT. These parameters, along with the associated standard^ 
deviations and F-ratios, (which in this case is (Parameter/a) ), 
are shown in Table 7.20. Simulations of these parameters for 
each of the optimal longitudinal stick input and the optimal 
collective stick input are shown in Figures 7. 22(a) -(f) and 7.23(a) 
respectively. These time histories are a marked improvement 
over the simulations of the initial parameters, although many 
discrepancies still exist. There are several possible reasons 
for this poor agreement. First, it is possible that there are 
unmodeled bias, scale factor and other measurement error effects 


in the 6 


COLL 


maneuver. Also, there may be substantial lateral 


coupling and unmodeled rotor effects present. Another possibil- 
ity is that key ^cqll Parameters had not converged resulting 

in the poor matches. This would be the case if the S r „ IT and 


the 


LONG 


COLL 

were at slightly different flight conditions, and 


the SCIDNT program converged to the 
of the higher signal level. 


LONG 


parameters because 


The SCIDNT results for Case 2 (doublet inputs) are presented 
in Table 7.21. Final parameter simulations to the longitudinal 
stick doublet input (Figure 7.17) and the collective stick doublet 
(Figure 7.18) are shown in Figures 7.24 and 7.25 respectively 
In this case, improved results were observed for the 

maneuver. Note that this does not necessarily mean that doublet 
inputs are "better" than optimum inputs. In order to properly 
compare doublet and optimal inputs, there must be a low level 
of external disturbances (i.e., gusts) and all other factors 
must be the same (particularly the energy of the input signals) . 

In addition, the pilot must be able to adequately duplicate the 
optimal stick input. Even then, the doublet inputs may give 
better results. The optimally designed stick input may not really 
be optimal because of the procedure used for designing that 
input. It is possible that modeling errors (coming about from 
higher degrees of freedom, nonlinearities and/or erroneous para- 
meter values) can cause large differences in the "optimal" inputs. 
In addition, there are inherent limitations in the multi-step 
input approach that can lead to differences. These limitations 
come about due to the size and number of steps in the designed 
input. Large steps (i.e., large AT ) can make identification 
of shorter period characteristics difficult. Increasing the 

170 



number of steps (to make the step size smaller) increases the 
computer time, core storage and accuracy requirements. Hence, 
there is also a limitation on the length of the maneuver. 


In order to obtain the best results at this flight condi- 
tion, an additional run was set up using the optimal ^lONG 


control input and the doublet 


6 COLL 


control input for simul- 


taneous processing. The results of this Case 3 are presented 
in Table 7.22 and the corresponding simulations in Figures 
7.26 and 7.27. These results are somewhat better but not 
completely converged. A comparison of the parameters estimated 
for all three cases are presented in Table 7.23. There are 
notable differences in these results which warrant further study 


One of the principal benefits of a complete state variable 
model is that many transfer functions can be rapidly obtained. 
For the Case 3 model identified from the flight data, the follow 
ing stick to pitch attitude transfer function is derived: 


8(s) = .17 s 2 +.100s+.0017 

6 LONG s 4 +1.18s 3 +1. 39s 2 +.10s+. 0044 


the poles of which are s^ ^ ~ ■•565 +_ j.99, s^ ^ = -.024 + j.182 

and the zeroes of which are -.560 and - .017. (Here ^lq^g is 

positive inches of aft stick.) The corresponding Bode plot is 
shown in Figure 7.28. 

Table 7.19 

Makeup of the Three Cases 



MANEUVER 

MANEUVER 


1 

2 

Case 1 

Optimal Longitudinal 

Optimal Collective 


Stick Input 

Stick Input 

Case 2 

Longitudinal Stick 

Collective Stick 


Doublet 

Doublet 

Case 3 

Optimal Longitudinal 

Collective Stick 


Stick Input 

Doublet 
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FIGURE 7.15(a).- OPTIMAL LONGITUDINAL STICK INPUT 

MANEUVER 
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FIGURE 7.15(b).- (CONTINUED) 
• S™ TT EXCITATION 
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FIGURE 7.16(a).- OPTIMAL COLLECTIVE 

STICK INPUT MANEUVER 
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FIGURE 7.16(b).- (CONTINUED) 
• 6„ nTT EXCITATION 
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FIGURE 7.17(a).- LONGITUDINAL STICK DOUBLET INPUT 




*-J / — > 

►J 2 

O t-i 
U'-' 
«o 


177 


oooc 



FIGURE 7.18(a).- COLLECTIVE STICK DOUBLET INPUT 

MANEUVER 



4.0 DUO 
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FIGURE 7.18(b).- (CONTINUED) 
* 6 COLL EXCITATI0N 
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INITIAL PARAMETER MODEL USED FOR UH-1H IDENTIFICATION 
(BASED ON C- 81 MODEL) 
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FIGURE 7.20(a).- SIMULATION OF INITIAL PARAMETER MODEL TO THE 

OPTIMAL LONGITUDINAL ST I OK INPUT COMPARED TO 
THE ACTUAL UH-1H RESPONSE 
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FIGURE 7.20(b).- (CONTINUED) 
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FIGURE 7.20(d).- (CONTINUED) 
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FIGURE 7.21(a).- SIMULATION OF INITIAL PARAMETER 

MODEL TO THE OPTIMAL COLLECTIVE 
STICK INPUT. COMPARED TO ACTUAL 
UH-1H RESPONSE 
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FIGURE 7.21(d).- (CONTINUED) 
o 0 MEASUREMENT 
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Table 7.20 

Parameter Estimates and Standard Deviations for UH-1H 
• Case 1 



INITIAL 

PARAMETER ESTIMATE 

FINAL PARAMETER 
ESTIMATE 

STANDARD 

DEVIATION 

F-RATIO 

x u , 1/sec 

-0.0227 

-0.0264 

0.00024 

12296.49 

Z u , 1/sec 

0.0013 

-0.0255 

0.00128 

394.00 

M u , 1/ft-sec 

0.0037 

0.0025 

0.00002 

11664.90 

X w , 1/sec 

0.0587 

0.0429 

0.00070 

3757.95 

Z w , 1/sec 

-0.7542 

-0.4258 

0.00312 

18572.47 

Mw, 1/ft-sec 

-0.0032 

-0.0093 

0.00005 

31510.22 

Mq, 1/sec 

-0.5305 

-0.7509 

0.00775 

9396.00 

x 6| nN R . ft/sec 2 - in 

0.8620 

0.5947 

0.0115 

2677.86 

Z «L0NG’ «/sec ! -in 

2.4600 

2.1250 

0.0602 

1246.41 

\0NG> 1/secM " 

-0.1800 

-0.1709 

0.0011 

25094.51 

x <Srm i » ft/sec 2 - in 

0.7115 

0.2051 

0.0094 

472.72 

z Srm i » ft/sec 2 - in 

-9.8180 

-2.9867 

0.0573 

2717.62 

m<s coll’ 1/sec2_in 

0.0064 

-0.0332 

0.0010 

1096.91 
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FIGURE 7.22(a).- SIMULATION OF CASE 1 FINAL PARAMETER 

MODEL TO THE OPTIMAL LONGITUDINAL 
STICK INPUT. COMPARED TO ACTUAL 
UH-1H RESPONSE 

* V TOTAL MEASUREMENT 
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FIGURE 7.23(a).- SIMULATION OF CASE 1 FINAL PARAMETER 

MODEL TO THE OPTIi'AL COLLECTIVE STICK 
INPUT. COMPARED TO ACTUAL UH-1H RESPONSE 
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FIGURE 7.23(b).- (CONTINUED) 
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FIGURE 7.23(d).- (CONTINUED) 
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FIGURE 7.23(e).- (CONTINUED) 
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Table 7.21 

Parameter Estimates and Standard Deviations for UH-1H 


• Case 2 



INITIAL 

PARAMETER ESTIMATE 

FINAL PARAMETER 
ESTIMATE 

STANDARD 

DEVIATION 

F-RATIO 

X u » 1/sec 

-0.0227 

-0.0328 

0.00063 

2713.21 

Z u , 1/sec 

0.0013 

-0.0221 

0.00478 

21.27 

M u , 1/ft-sec 

0.0037 

0.0040 

0.00007 

3430.32 

X w , 1/sec 

0.0587 

0.0257 

0.00201 

162.97 

Z w , 1/sec 

-0.7542 

-0.5268 

0.00601 

7681.96 

Mw, 1/ft-sec 

-0.0032 

-0.0075 

0.00008 

8471.64 

M q , 1/sec 

-0.5305 

-0.5170 

0.01405 

1354.68 

VoNG 1 ft/S “ 2 - in 

0.8620 

0.4725 

0.01556 

921.96 

Z «L0NG' ft/SecM " 

2.4600 

2.7983 

0.09074 

951.04 

VnG 1 ' /Sec2 - in 

-0.1800 

-0.1103 

0.00151 

5356.71 

x 6 C oll* ft / sec2 ~ in 

0.7115 

0.2907 

0.01036 

788.00 

ZficoLL* ft/ sec 2 - in 

-9.8180 

-5.1817 

0.07954 

4243.93 

M6 coll’ ,/sec! - in 

0.0064 

-0.0583 

0.00130 

2001.03 
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FIGURE 7. 24(a) .-SIMULATION OF CASE 2 FINAL PARAMETER 

MODEL TO THE LONGITUDINAL STICK DOUBLET 
INPUT. COMPARED TO ACTUAL UH-1H RESPONSE 

* V TOTAL MEASUREMENT 



Measured 

Simulated 







Measured 

Simulated 



209 


FIGURE 7.24(c).- (CONTINUED) 
• q MEASUREMENT 
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I FIGURE 7.25(a).- SIMULATION OF CASE 2 FINAL PARAMETER 

MODEL TO THE COLLECTIVE STICK DOUBLET 
INPUT. COMPARED TO ACTUAL UH-1H RESPONSE 
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FIGURE 7.25(c).- (CONTINUED) 
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TABLE 7.22.- PARAMETER ESTIMATES AND STANDARD DEVIATIONS FOR UH-1H 


• Case 3 



INITIAL 

PARAMETER 

ESTIMATES 

FINAL 

PARAMETER 

ESTIMATES 

STANDARD 

DEVIATION 

F-RATIO 

X u , 1/sec 

-0.0227 

-0.0271 

0.00026 

11092.39 

Z u , 1/sec 

0.0013 

-0.0226 

0.00144 

244.79 

H , 1/ft-sec 

0.0037 

0.0027 

0.00003 

9694.58 

x w , 1/sec 

0.0587 

0.0391 

0.00075 

2701.15 

Z w » 1/sec 

-0.7542 

-0.4573 

0.00337 

18361.13 

M w , 1/ft-sec 

-0.0032 

-0.0095 

0.00006 

28760.19 

M q , 1/sec 

-0.5305 

-0.7017 

0.00697 

10141.01 

X 6L0NG’ ft/sec 2 - in 

0.8620 

0.5871 

0.01119 

2754.05 

Z 6L0NG’ f t/sec 2 -in 

2.4600 

1.9616 

0.06012 

1064.63 

M 6L0NG' ’/secMn 

-0.1800 

-0.1685 

0.00106 

25282.50 

X 6C0LL’ ft/sec2 - ,n 

0.715 

0.1691 

0.01151 

215.80 

Z 6C0LL- ft/secM " 

-9.8180 

-3.9439 

0.08663 

2072.84 

M 6C0LL- '/secMn 

0.0064 

-0.0535 

0.00163 

1072.94 
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FIGURE 7.26(a).- SIMULATION OF CASE 3 FINAL PARAMETER 

MODEL TO THE OPTIMAL LONGITUDINAL STICK 
INPUT. COMPARED TO ACTUAL UH-1H RESPONSE 

• MEASUREMENT 
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FIGURE 7.26(b).- (CONTINUED) 
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FIGURE 7.26(e).- (CONTINUED) 
• a MEASUREMENT 
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FIGURE 7.27(a).- SIMULATION OF CASE 3 FINAL PARAMETER 

MODEL TO THE COLLECTIVE STICK DOUBLET 
INPUT. COMPARED TO ACTUAL UH-1H RESPONSE 

* V TOTAL MEASUREMENT 
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FIGURE 7.27(b).- (CONTINUED) 
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TABLE 7.23.- UH-1H FLIGHT IDENTIFIED PARAMETERS - COMPARISON OF THE THREE CASES 


• 60 Knots 

• 2 Maneuver Runs 



INITIAL 

PARAMETER 

ESTIMATES 

CASE 1 

FINAL PARAMETER 
ESTIMATES 

CASE 2 

FINAL PARAMETER 
ESTIMATES 

CASE 3 

FINAL PARAMETER 
ESTIMATES 

X u , 1/sec 

-0.0227 

-0.0264 

-0.0328 

-0.0271 

Z u , 1/sec 

0.0013 

-0.0255 

-0.0221 

-0.0226 

M u> 1/ft-sec 

0.0037 

0.0025 

0.0040 

0.0027 

X w , 1/sec 

0.0587 

0.0429 

0.0257 

0.0391 

Z w , 1/sec 

-0.7542 

-0.4258 

-0.5268 

-0.4573 

M w , 1/ft-sec 

-0.0032 

-0.0093 

-0.0075 

-0.0095 

M q , 1/sec 

-0.5305 

-0.7509 

-0.5170 

-0.7017 

x ai nw r» ft/sec 2 -in 
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Z 6C0LL* ft/ 5 ** 2 - 1 " 
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M SC0LL’ '/sec 2 -in 
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CHAPTER VIII 
CONCLUSIONS 


An integrated methodology for rotorcraft system identifi- 
cation has been described. This methodology consists of three 
distinct data processing steps and a technique for designing 
inputs to improve the identif iability of the data. These 
elements are as follows: 

(1) A Kalman Filter/Smoother Algorithm which 
estimates states and sensor errors from 
error- corrupted data. Gust time histories 
and statistics may also be estimated. 

(2) A Model Structure Estimation Algorithm for 
isolating a model which adequately explains 
the data. 

(3) A Maximum Likelihood Algorithm for estimating 
the parameters and estimates for the variance 
of these estimates. 

(4) An Input Design Algorithm, based on a maximum 
likelihood approach, which provides inputs 

to improve the accuracy of parameter estimates. 

A discussion of each step is presented, with examples to 
both flight and simulated data cases. 

Simulated data examples indicate that the software is 
valid for the idealized errors of such synthetic tests. 

Such simulation data processing provides intuitive insight 
into the manner in which various errors affect the identifi- 
cation performance. One unexpected benefit of applying the 
algorithms presented as part of the methodology is that 
simulation bugs were rapidly found and, in some cases, 
corrected . 

The following principal conclusions are made from the 
experience gained on the flight data processing. 

(1) Any set of flight data poses unique problems 
for estimation of the data errors, model 
structure estimation, and parameter identifi- 
cation. Maximum flexibility on the part of 
software is needed to quantify the numerous 
possibilities of errors. 

(2) Rotorcraft can pose a particularly severe 
instrumentation problem. For the data 


235 



provided to us for this effort, no rotorcraft, 
as expected, had a complete set of perfect 
sensors. Sensors were either completely 
absent, failed, or badly degraded on the 
data. This stresses the importance of the 
state estimation aspect of the methodology, 
which can be used to reconstruct very poor 
data channels. 

Use of a sequential processing scheme provides 
valuable engineering judgement to be applied 
at key points in the processing. The three 
basic software algorithms can be used to 
actually check and double-check some estimates 
of each other. For example, sensor bias can 
be estimated by all the algorithms independently. 
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APPENDIX A 

MAXIMUM LIKELIHOOD WITH NO PROCESS OR MEASUREMENT NOISE 


The maximum likelihood method can be simplified when 
either process noise or measurement noise are absent. 

No process noise .— If the process noise is zero and 
initial states are known perfectly, i.e., w(t) and P(0) are 
zero, the covariance of the error in the predicted state is 
also zero. It is clear that Kalman gains are zero. The 
innovations are the output error, i.e., 

v(i) = y (i) - h (x(t i ) ,u(t i ) ,0,t i ) (1) 

and the innovation covariance is 

B(i) = R (2) 

the log-likelihood functionis. 


Log (i?Ce|z)) “-is v (i) R~ v (i) + log |R| (3) 

2 i=l 

which on optimizing for unknown parameters in R gives 

R = i S v(i) v T (i) (4) 

N i=l 

The equality in (4) holds only for those elements of R which 
are not known a priori. For instance, even if R is known to 
be diagonal, the right hand side matrix will not be diagonal, 
in general; but, the off-diagonal terms should be ignored 

A 

before they are equated to R. Using (4) in (3) 

N 

Log (^(8 |z)) = - j E v^(i) R v(i)n+ constant (5) 

2 i=l 

The optimizing function is the same as that for the output 
error method except that the measurement noise covariance 
matrix is determined using (4) and is used as the weighting 
matrix in the criterion function. In the output error method, 
the measurement noise is assumed known and the weighting 
function is arbitrary. 
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The first and second derivatives of the log-likelihood 
function with respect to unknown parameters are 


3 

30- l0g 

(£?(Q 1 z) ) = 

N 

- Z 
i=l 

v T Ci) R' 1 3 ^> 
39 3 

(6) 

3 2 log (^(0 lz) ) _ 

N , 
y i 

\ 3v T (i) --1 3v ( i) 


30 . 
J 

a? 

a> 

l 

i-i ! 

1 3e k 39 i 



* v T (i) 5-1 j 

J 30, 30, [ 


3 k 

The terms in the second derivative are approximated as 

.2 

-tog 

^k 


3 log (jgce lz)) _ 
30 j 30 , 


N 

Z 

i=l 


3v T (i) n-1 3v(i) 
30 k 30 j 


(7) 


( 8 ) 


No measurement noise .— If all states are measured with 
no noise, the covariance of the error in state estimates is 
zero at the beginning of any time update. 


P(i-1 1 i-1) = 0 
and x(i-l|i-l) = x(i-l) 


( 9 ) 


It is easy to show in this case that for fast sampling the 
log-likelihood ' function is quadratic in the difference be- 
tween measured values of x and f(x,u,0,t). The method re- 
duces to the equation error method, the weight W being chosen 
as 


W = 


1 

T 



f (x,u,0 ,t)) 



f (x,u 



( 10 ) 


Thus, the maximum likelihood method and equation- error 
methods are equivalent except for the technique for choosing 
the weighting matrix. 
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APPENDIX B 


This Appendix gives a brief description of Walsh functions and 
a short summary of their properties. The Walsh functions <j> Q (t), 

(t) ,<{>•• • are a set of square waves which are orthonormal. 

Each Walsh function can be decomposed into more elementary square 
waves, or Rademacher functions. 

Rademacher functions r. (t) , are a set of square waves of 

K fl-kl 

unit height with periods equal to 1,1/2, 1/4, . . . , 2 : J In other 

words, the number of cycles of the squre waves of r v (t) is 

k-1 K 

2. A few Rademacher functions are shown in Fig. B.l. The Walsh 

functions are defined in terms of r^ft) as follows: 

♦o (t) ■ r 0 Ct) 

^(t) = r x (t) 

9 

^k ^k-1 ^k-2 

*i Ct) = [r k (t ^ ] * f r k-l^^ 

where 

k = [log 2 i ] + 1. 

[•] means taking the integral part and • • • * b^ is 

the binary number expression of i. Typical Walsh functions are 
shown in Fig . B . 2 . 

Next, we consider some of the properties of the Walsh func- 
tions . 

Integration of Walsh Functions .— Integrals of Walsh functions 
may be approximated by sums of Walsh functions. It is shown in 
Reference [4] that for the first four functions 




Higher dimension P matrices may be obtained straightforwardly. 


Evaluation of <j>^ (t) $ (t) If two Walsh functions are multi- 

plied together the produce is a Walsh function obtained by the mod 
2 addition of the binary representation of the original functions. 
Therefore, <j>^(t)$(t) may be written in terms of Walsh functions. 

It is easy to show that 


4> 0 (t) $ (t) = 


0 

1 

0 

0 


0 

0 

1 

0 


Kt) 


<j> 1 (t)$(t) 


0 

1 
0 
0 


1 

0 

0 

0 


0 

0 

0 

1 


0 

0 

1 

0 


<Kt) 


(t) = 


0 

0 

0 

1 


1 

0 

0 

0 


0 

1 

0 

0 


*(t) 


4>3 (t)*(t) = 


0 

0 

0 

1 


0 

0 

1 

0 


0 

1 

0 

0 


1 

0 

0 

0 


*(t) 


Delay Matrix.— A delayed 4>(t) may also be written in terms 
of the Walsh functions. It is easy to verify that for first four 
Walsh functions 


Ct-r) 



*00 


- D 4x4*M 
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It is clear that 


will produce two units of delay, e.g.. 




= D $(t) = 


1 

16 


-8 

-8 

0 

0 


0 0 

0 0 

8 -8 

8 -8 


*(t) 
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APPENDIX C 


In this appendix we derive expressions for the information 
matrix and the quadratic constraint when the state and outputs are 
assumed to be multistep. Walsh functions are used in this deriva- 
tion. Details of the properties of the Walsh functions are cov- 
ered in References [56], [57], and [58]. The important properties 

for the purpose of the present application are summarized in 
Appendix . To simplify the various derivations it is assumed 

that s = 2^, where Q is an integer. Let cf>^, i=0 , 1 , . . . , s - 1 be 

the ith Walsh function in the dyadic order and let the input be 
assumed to be 


u(t) = H (t) (1] 

$(t) is a vector of Walsh functions from 0 to s-1. Also, let 
* CO - W $ (t) (2) 

Since x(0)=0, 

x(t) = W f $(t) dr = W P $ (t) (3) 

J 0 

where P is an sxs matrix described^in Appendix A. Next, we 
approximate time varying matrices A, B, and C by Walsh series 

s-1 

A (t) = Z A C i ) <0(0 
i=0 

s-1 

B (t) = Z B (i) <j> . C t ) (4) 

i = 0 1 

s-1 

T k (t) = Z T k (i) <t» i (t) 


Using this approximation. Equation (6) becomes 


s-1 

W$(t) = z A(l)c() . (t) WP$(t) 
i = 0 x 


s-1 

+ l BCi)<j>, CO H«CO (5) 
i = 0 


where U(i) is an sxs matrix such that (see Appendix A) 
4> i (t)®(t) = U(i)$(t) 
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Matrix H can be rearranged into a vector as before and the matrix 
W may be rearranged into another n(m+l)s vector w by substi- 
tuting its first column into the first n(m+l) element of w, the 
second column into the next n(m+l) elements and so on. Then 
Equation (5) becomes 

| I - V A(i) ©(PU(i)) T > w = V B (i)©U T (i) h (6) 

( i-0 ) i=0 

where 


A © B = 


therefore 


B 11 A 


B 21 A 


B 12 A 


B 22 A 


B .A B 0 A 
nl n2 


B, A 
In 


B 2n A 


B A 
nn 


(7) 


s-1 

W = I -{ E A(i) © (PU(i)) 
(i=0 

A L h 


1-1 s-1 T 

’ Z B(i) © U 1 (i) h 
1 i=0 


( 8 ) 


Similarly, Equation (6.18) may be written as 
9 z s-1 

96, = E T,(i) <j>.(t) W*(t) 

K i=0 K 1 

(9) 

s-1 

= Z T, (i) WP U(i)*(t) 
i = 0 K 


A typical element of the information matrix is 

s-1 


M 


T , . 


T t .,T 


kS. 


i, j=l 


Tr T/(i) R x T 9 ( j ) WP U(j) IT (i) P^W 1 (10) 


Since w is a linear function of h, M k£ may be written as 


M k£ = h* V(k , 2.) h 


(ID 
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An explicit expression for V (k , £) may be obtained by the reader. 
The quadratic constraint may also be written in terms of the vector 
h. 


s -1 

£ T n (i) WPU(i) $ (t) 
i=0 u 


T 0 (j)WPU(j)*(t) 


+ u T (t)u(t) jdt 

I s - x 

jr .2 Tq (i)Q T q C j ) WPU(j) U T (i) P T W T 
+ H T H 1 = h T Q'h (12) 


We note that in both cases (see 6.2.2) similar kinds of relation- 
ships are obtained for the information matrix and the state and 
input quadratic constraint. 
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